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I.   INTRODUCTION 

A.       GENERAL 

Any  study  of  the  motion  of  a  ship  through  water  requires  the  use  of  models. 
The  model  could  be  a  mathematical  representation  which  attempts  to  describe  the 
situation  symbolically,  searching  for  a  quantitative  solution,  or  a  physical  model 
which  strives  to  replicate  or  simulate  the  actual  conditions  for  the  purpose  of 
collecting  empirical  data.  The  two  can  be  combined  such  that  they  complement  each 
other,  or  most  often,  physical  models  are  used  to  test  the  validity  and  accuracy  of 
mathematical  models.  However,  both  models  fail  to  completely  describe  the  physical 
phenomenon.  The  physical  model  falls  short  by  the  inability  to  simultaneously  satisfy 
all  the  parameters  describing  the  problem.  The  mathematical  model  suffers  by 
departures  from  the  true  state  through  assumptions  made  during  the  development 
of  the  equations  [Ref.  1].  Notwithstanding  these  shortcomings,  the  use  of 
models  in  predicting  a  ships  response  to  its  environment,  i.e.,  wind  and  waves,  and 
its  own  propulsive  and  control  systems  is  essential  in  the  design  process  to  ensure 
safety  of  the  ship,  performance  of  its  mission,  its  survivability  in  extreme  conditions 
and  its  efficiency  during  normal  operating  conditions.  [Ref.  1] 

In  general,  the  examination  of  ship  motions  is  separated  into:  1)  steering  and 
maneuverability  (calm  water),  and  2)  seakeeping  (waves,  current),  with  motion 


stability  (absence  of  external  excitation)  and  control  (external  excitation)  some  of  the 
primary  concerns.  [Ref.  1] 

The  basis  of  the  mathematical  model  is  rigid  body  dynamics,  hence  Newton's 
laws  of  motion:  [Ref.  2] 

d 


F    =    —  x  {linear  momentum) 

*     h  (ii) 

M    =    —  x  (angular  momentum) 
dt 


The  ship  is  free  to  rotate  (roll,  pitch,  yaw)  and  translate  (surge,  sway,  heave)  in  all 
six  degrees  of  freedom  and  the  forces  and  moments  acting  on  the  vehicle  are 
comprised  of:  [Ref.  2] 

1.  Hydrodynamic  forces  and  moments  on  the  bare  hull,  appendages,  rudder  and 
propeller. 

2.  Inertial  reaction  forces  and  moments. 

3.  Wind,  waves  and  currents. 

4.  External  forces. 

Moreover,  these  forces  and  moments  are  dependent  upon  properties  of  the  rigid 
body  (e.g.,  geometry,  mass,  center  of  gravity),  its  orientation  (i.e.,  to  inertial  frame 
of  reference  and  to  the  body  fixed  frame  of  reference),  its  dynamics  (velocities, 
accelerations,  propeller  speed,  rudder  deflection  and  deflection  rate),  and  the  fluid 
properties  (e.g.,  density,  viscosity,  pressure,  energy).  [Ref.  1] 


At  this  point,  the  mathematical  model  has  rapidly  developed  into  a  very 
complex  and  cumbersome  system  and  requires  simplification  in  order  to  be  of 
practical  use.  Choice  of  the  body  fixed  coordinate  system  to  take  advantage  of  any 
symmetries  is  one  simplification  used  almost  exclusively.  Eda,  et  al,  [Ref.  3] 
found  in  one  study  of  bulbous  bow  type  ships  that  the  slight  asymmetry  of  the 
underwater  hull  form  increased  significantly  in  roll  such  that  this  simplification  was 
not  valid.  Another  reduction  in  complexity  often  used  is  to  separate  the  lateral 
motions  (sway,  yaw,  roll)  from  the  longitudinal  motions  (surge,  pitch,  heave) 
[Refs.  4,  5,  6].  Later  studies  [Refs.  3,  7,  8,  9,  10]  account  for  the  loss 
of  speed  during  a  turn  by  including  the  surge  equation  with  the  lateral  motions  for 
maneuvering  studies.  Seakeeping  studies  [Refs.  8,  11,  12,  13]  heretofore  have 
separated  out  the  roll  equations  using  a  one  degree  of  freedom  system  only,  coupling 
effects  were  neglected  as  very  small  or  incorporated  as  forcing  terms.  Rutgersson  and 
Ottosson  [Ref.  8]  superimposed  the  maneuvering  model  motions  with  the  seakeeping 
motions. 

As  can  be  ascertained  from  this  discussion,  the  forces,  moments,  velocities  and 
accelerations  acting  upon  the  vessel  are  highly  nonlinear,  involving  complex  coupling 
of  terms.  Linear  theory  has  been  successful  in  predicting  and  analyzing  directionally 
stable  ships  with  controls  fixed  for  small  perturbations  [Ref.  14].  In  addition,  the 
linear  hydrodynamic  forces  and  moments  have  been  reduced  to  semi-empirical 
equations  based  on  ship  design  characteristics  such  as  block  coefficient,  length, 


breadth  and  draft  [Refs.  2,  15].  However,  there  has  been  no  completely  analytical 
procedure  which  predicts  the  nonlinear  hydrodynamic  forces  and  moments  [Ref.  2]. 

It  is  necessary  therefore,  to  combine  the  mathematical  model  with  the  physical 
model  such  that  the  physical  model  is  used  to  measure  these  nonlinear  forces  and 
moments  and  incorporate  the  results  in  the  mathematical  model.  "The  functions  that 
may  approximate  the  hydrodynamic  forces  and  moments  can  be  expressed  formally 
by  Taylor  expansion  around  the  state  of  equilibrium  with  respect  to  the  quantities 
affecting  the  different  motions,  such  as  the  axial  speed  u,  turning  rate  r,  rudder  angle 
8,  etc.  The  polynomial  coefficients  thus  describing  the  hydrodynamic  forces  and 
motions  can  be  determined  from  captive  model  tests  at  different  kinds  of  basic 
motions  [Ref.  8]. "  A  complete  treatment  of  the  captive  model  test  procedure  can  be 
found  in  Principals  of  Naval  Architecture,  Vol.  Ill  [Ref.  2],  and  will  not  be  discussed 
here.  The  Taylor  expansion  and  subsequent  model  development  is  elaborated  in  the 
following  Chapters. 

Captive  model  tests  for  the  lateral  motions  are  more  difficult  than  for  the 
longitudinal  motions  due  to  the  large  sizes  required  to  adequately  represent  the 
physical  model.  Consequently,  there  exists  little  data  for  the  third  order  coefficients 
necessary  to  produce  semi-empirical  equations  for  the  nonlinear  hydrodynamic  forces 
and  moments  arising  from  these  lateral  motions  similar  to  those  which  exist  for  the 
first  order  terms.  One  study  by  Inoue,  et  al,  [Ref.  15]  has  produced  formulations  for 
limited  nonlinear  terms  in  yaw  and  sway  based  on  drift  angle  and  turning  rate.  The 
yaw  moment  coefficients  so  produced  are  not  in  good  agreement  with  the  model  test 


data  under  all  conditions  of  loading  and  motions  induced.  Another  study  done  by 
Son  and  Nomoto  [Ref.  9],  utilizes  model  test  data  for  the  third  order  hydrodynamic 
coefficients  of  a  SR-108  container  ship.  In  the  absence  of  such  semi-empirical 
formulations  for  the  forces  and  moments  in  the  equations  of  motion  and  access  to 
model  test  facilities,  the  aforementioned  data,  together  with  vehicle  design 
characteristics,  can  be  used  in  a  mathematical  model  to  ascertain  particular 
information  about  the  coupling  of  yaw,  sway  and  roll  motions. 

B.      SCOPE  OF  THIS  THESIS 

This  research  is  intended  as  a  bridge  between  the  static  roll  restoring  moment, 
controls  fixed  analysis  of  the  coupled  lateral  motion  problem  and  a  fully  coupled 
treatment  of  a  ships  roll  response  to  maneuvering  in  a  seaway.  As  discussed 
previously,  current  studies  in  the  area  of  maneuvering  and  steering  range  from  three 
degree  of  freedom  models,  to  four  degree  of  freedom  models  which  account  for  the 
speed  loss  during  a  turn.  These  studies,  however,  use  the  equilibrium  point  in  the 
Taylor  expansion  of  the  representative  hydrodynamic  forces  and  moments,  as  straight 
upright,  zero  rudder  deflection,  such  that  v0=r0=<£0=50=0  and  all  perturbations  and 
nonlinearity  effects  are  evaluated  based  on  this  single  nominal  point.  This  research 
has  as  its  basic  premise  that  the  nominal  point  for  the  Taylor  expansion  changes 
significantly  with  rudder  angle.  As  a  result  of  this  change,  the  stability  characteristics 
of  the  ship  may  also  change  appreciably,  a  fact  which  needs  to  be  established  for 
successful  prediction  of  ship  maneuvering  response  in  a  seaway. 


Therefore,  the  objective  is  to  determine  the  significance  of  the  coupling  of  roll 
into  sway  and  yaw  motions  as  a  function  of  rudder  angle.  In  analyzing  this  coupled 
problem,  the  effects  of  roll  into  sway/yaw  and  conversely,  sway/yaw  into  roll  were 
investigated. 

Development  of  the  mathematical  model  for  the  coupled  yaw,  sway,  and  roll 
equations  of  motion  along  with  models  for  the  rudder  forces  and  moments,  the 
righting  moment,  and  the  hydrodynamic  forces  and  moments  is  presented  in  Chapter 
II.  The  final  form  of  the  model  being  reduced  to  a  set  of  four  nonlinear,  coupled, 
first  order  differential  equations.  The  additional  equation  is  derived  from  kinematics 
to  relate  the  roll  angle  and  the  roll  rate. 

In  order  to  utilize  Taylor  expansion  for  formulation  of  the  functions  which 
approximate  the  hydrodynamic  forces  and  moments,  the  state  of  equilibrium  for  the 
motion  parameters,  v,  r,  and  <f>,  is  established.  This  is  accomplished  through  steady 
state  analysis.  In  Chapter  III,  this  analysis  is  first  done  for  the  linear  case  to  obtain 
an  initial  approximation  for  the  case  of  near  zero  rudder  deflection  (8—0).  This 
linear  steady  state  solution  is  incorporated  into  the  nonlinear  steady  state  analysis 
in  an  iterative  process  which  produced  the  equilibrium  values  for  yaw  rate,  sway 
velocity,  and  roll  angle  as  a  function  of  the  rudder  deflection  in  a  steady  turn.  As  the 
rudder  angle  is  incrementally  increased,  the  steady  state  values  from  the  previous 
rudder  angle  are  used  as  the  initial  approximations  to  initiate  the  iterations  for  the 
current  rudder  deflection. 


Once  the  steady  state  solutions  have  been  computed,  their  stability  properties 
are  established  by  local  perturbation  of  the  coupled  equations  of  motion  in  the 
neighborhood  of  each  nominal  point.  This  procedure,  described  in  Chapter  IV,  yields 
a  generalized  eigenvalue  problem.  The  solutions  to  the  eigenvalue  problem  provide 
information  as  to  the  stability  properties  of  the  nominal  points  and  thus  characterizes 
the  stability  of  the  coupled  sway,  yaw  and  roll  motions.  Demonstration  of  the 
coupling  effect  on  the  lateral  motions,  is  shown  utilizing  comparison  to  two 
independent  systems: 

1.  A  sway/yaw  model  represented  by  two  nonlinear  equations  decoupled  from 
roll  in  that  the  hydrodynamic  forces  and  moments  are  functions  of  sway 
velocity  (v)  and  yaw  rate  (r)  only. 

2.  A  roll  model  represented  by  the  nonlinear  roll  equation  and  the  kinematic 
roll  rate  equation  decoupled  from  sway  and  yaw  in  that  the  hydrodynamic 
forces  and  moments  are  functions  of  roll  angle  (<£)  and  roll  rate  (p)  only. 

The  directional  and  roll  stability  characteristics  in  a  steady  turn  are  established 
by  the  degree  to  which  the  real  part  of  the  complex  eigenvalues  is  negative  and  the 
magnitude  of  the  roll  damping  ratio. 

Chapter  V  relates  the  sensitivity  and  stability  of  design  parameters  such  as 
metacentric  height  and  trim  to  the  coupling  between  the  lateral  motions.  Results  are 
presented  in  terms  of  rudder  angle  with  a  qualitative  summary  at  the  end  of  this 
chapter. 

Finally,  Chapter  VI  summarizes  the  conclusions  and  presents  recommendations 
for  further  research  and  utilization  of  the  mathematical  model  herein  developed. 


II.   MATHEMATICAL  MODEL 

A.       EQUATIONS  OF  MOTION 

For  ship  motions,  the  full  nonlinear  equations  of  motion  in  a  body  fixed  axis 
system  are  shown  in  Equations  (2.1).  These  equations  reflect  an  assumption  of 
symmetry  along  the  longitudinal  axis  (i.e.,  yG=0).  The  orthogonal,  right-hand  axes  for 
the  body  fixed  system  is  shown  in  Figure  2.1.  The  parameters  are  defined  in  Table 
I. 

Surge:      m[u-x(J(q2+r2)+z(J(.q+pr)-rv+wq]  =  X  (2.1a) 

Sway:       m[v+xG(pq+f)+zG(qr-p)+ur-wp]  =  Y  (2.1b) 

Heave:      mlw+x^rp-fy-ZQ^+q^+pv-qu]  =  Z  (2.1c) 

Roll:  z    ^       *>  *  (2.1d) 

=  K-AGZ($) 

Pitch:  lzie; 

Xgiw+pv-qu)]  =M 

Yaw:       If+(L-Ix)pq+Ixz(qr-p)+mxG(v+ur-wp)  =  N  (2.  If) 


Figure  2.1:   Coordinate  System  [Ref.  11] 
B.       SIMPLIFYING  ASSUMPTIONS 

For  the  purpose  of  this  investigation,  a  mathematical  model  having  three 

degrees  of  freedom  vice  the  six  degrees  of  freedom  depicted  in  Equations  (2.1)  was 

used.  The  following  simplifying  assumptions  were  made: 


1.  The  rotational  velocity  and  acceleration  about  the  y-axis  are  zero.     (  q  =  0 
and   q  =  0  ) 

2.  The  translational  velocity  and  acceleration  in  the  z  direction  are  zero.     (  w 
=  0   and   w  =  0  ) 

3.  The  vertical  heave  and  pitch  motions  are  decoupled  from  the  horizontal  plane 
motions. 

4.  The  product  of  inertia  1^  is  very  small  and  can  be  neglected. 

5.  The  surge  equation  is  substituted  by  an  algebraic  equation  which  is  a  function 
of  u,  V,  and  5. 

6.  The  longitudinal  center  of  gravity,  (LCG)  and  the  longitudinal  center  of 
buoyancy,  (LCB)  are  at  midship. 


7.  The  vertical  center  of  gravity,  (VCG)  is  on  the  centerline. 

8.  The  only  important  forces  and  moments  acting  on  the  ship  induced  by  the 
rudder  are  those  due  to  rudder  deflection.  Forces  and  moments  due  to  8  and 
8  are  negligible.  [Ref.  2] 

TABLE  I.   EQUATIONS  OF  MOTION  PARAMETERS 


Parameter 

Description 

x,y,z 

Distance  along  the  principal  axes 

u,v,w 

Translational  velocity  components  of  ship  relative  to  fluid 
along  body  axes 

p>q>r 

Rotational  velocity  components  of  ship  relative  to  inertial 
reference  system  along  body  axes 

X,Y,Z 

Hydrodynamic  force  components  along  body  axes 

K,M,N 

Hydrodynamic  moment  components  along  body  axes 

ip,6,<f> 

xp  yaw  angle:  bow  to  starboard  positive 

0  pitch  angle:  bow  up  positive 

<f>  roll  angle:     starboard  down  positive 

m 

Mass  of  ship 

XgiYg'^g 

Coordinates  of  the  center  of  gravity  in  the  body  axis 
system 

Ix>Iy>Iz 

Moments  of  inertia  about  the  body  axis  system 

Ixz>Iyz>Ixy 

Products  of  inertia  about  the  body  axis  system 

V 

Displacement  volume  of  ship 

A 

Displacement  weight  of  ship 

GZ(4>) 

Righting  moment  as  a  function  of  roll  angle 

8 

Rudder  angle  in  radians 

V 

Initial  velocity  of  ship 

P 

Mass  density  of  sea  water 

L 

Ship  length  between  perpendiculars  (LBP) 

10 


Applying  these  assumptions  to  Equations  (2.1),  the  three  degree  of  freedom 
equations  of  motion  for  the  model  used  are  given  in  Equations  (2.2). 


Sway:         m[v+x(/-z(^+ur]  =  Y 


(2.2a) 


Roll:      /^-mz^v+wr)  =  tf-AGZ(4>) 


Yaw:      Ij+mxjiy+ur)  =  N 


(2.2b) 
(2.2c) 


Equations  (2.2)  are  nondimensionalized  for  ease  of  working  between  model  test 
data  and  actual  ship  test  data  using  the  relationships  shown  in  Table  II. 

To  demonstrate  the  use  of  these  nondimensionalizing  terms,  substitution  of  the 
appropriate  values  from  Table  II  into  Equation  (2.2a)  yields: 


(1 


**• 


m'[ 


(      o\ 

V2 


U 


v'+(Lxc) 


(v2  \    ,iv  ;i 


^(rMi"1^) 
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TABLE  II.   NONDIMENSIONAL  PARAMETER  RELATIONSHIPS 


,'    - 


I    r  \ 


V1      =      VI 


\Y    ) 


Y'    =    Y 


12K         , 


■  ■  <*: 


j   _ 


r'    =    r 


<L* 


V2) 


N'    =    N 


2P 


♦'    =    <l> 


*'    = 


a 


-I   _ 


K'    =    K 


'       1 


±pL*V2 


P      =    P- 


[V) 


P'    =    P 


,/  _ 


u 


m'    =    m 


<  i  N 


V 


12 


xa    =    * 


® 


zG    =    z 


® 


',    •    ', 


'  1  * 


>5, 


'*  -  A 


'  1  N 


dpL5 


AGZ(4>)'  -- 
AGZ((|)) 


'       1       N 


IpL3^2 
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simplifying  by  factoring  out  (V2/!.): 

f-pL2vA  m'  [v'+x'j'-z&'+u'r']    =    (-pL2vA  Y' 

and  Equation  (2.2a)  becomes  Equation  (2.3)  in  nondimensional  form. 

m'[v'  +  x'Gr'  -z'Gp'  +  u'r' ]  =  Y'  (2-3) 

Since  Equation  (2.3)  and  Equation  (2.2a)  are  of  the  same  form,  the  prime  notation 
will  be  dropped  and  all  equations  will  be  considered  as  represented  in 
nondimensional  form  unless  otherwise  indicated. 

C.      FORCE  AND  MOMENT  REPRESENTATION 

Using  the  method  of  Abkowitz  and  Strom-Tejsen  [Ref.  2],  the  sway  force,  roll 
moment  and  yaw  moment  can  be  expressed  as  Equations  (2.4). 

Sway  Force:         Y    =   /j(K,v,r,w,v,r,<|>,(i>,6)  (2-^a) 

Roll  Moment.       K    =   /2(«,v,r,w,v,r,<t>,<i>,$,6)  (2.4b) 

Yaw  Moment:         N    =   /3(K,v,r,ii,v,r,<t>,<j>,5)  (24c) 

A  third  order  Taylor  expansion  of  fv  f2,  and  f3  can  be  expressed  as  Equations 
(2.5). 

f,  =  Yvv+Y,v+Y^+Yvrrvr2+Y/+Yrnr*+Y„v2r+ 

(2.5a) 

Y^+Y^+Y^v2<b+Y^2+Y^rtf+Y^r24>+Y(6)+Y0 
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^v2<l>+^W)2^^r<t>2^(6)^0-A  GZ(4>) 

(2.5c) 

The  coefficients  in  Equations  (2.5)  are  the  hydrodynamic  coefficients  and 
represent  the  partial  derivative  with  respect  to  the  subscripted  variable.  For  example, 

(&y\ 

¥-L  .  The  terms  Y0,  Ko,  and  N0  are  the  sway  force, 
{dr>) 
roll  moment  and  yaw  moment  induced  by  the  propeller  respectively.  In  a  similar 


Yv  means  —  and  YvrT  means  — 
dv  dv 


manner  Y(5),  K(S)  and  N(5)  represent  the  force  and  moments  induced  by  the 
rudder.  The  righting  moment  due  to  the  static  stability  of  the  ship  is  represented  by 
the  AGZ(<£)  term. 

D.      RUDDER  FORCE  AND  MOMENT  REPRESENTATION 

The  expressions  used  to  determine  the  rudder  force  and  moments  were  taken 
from  Son  and  Nomoto  [Ref.  9]  and  are  presented  in  Equations  (2.6)  and  (2.7).  The 
parameters  are  defined  in  Table  III. 

r(6)  =    -  (1  +  aH)FNcos(6)  (2.6a) 

K(b)    =       (1  +  aH)zRFNcos(b)  (2.6b) 

ATO    =    "  (xR  +  ^tf)FNcos(5)  (2.6c) 
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TABLE  III.   RUDDER  FORCE  AND  MOMENT  PARAMETERS 


Parameter 

Description 

aH 

Rudder  to  hull  interaction  coefficient 

FN 

Normal  force  action  on  the  rudder 

Zr 

z  coordinate  of  point  on  which  rudder  force  Ys  acts 

Xr 

x  coordinate  of  point  on  which  rudder  force  Y5  acts 

xh 

x  coordinate  of  point  on  which  normal  force  FN  acts 

A 

Rudder  aspect  ratio 

AR 

Rudder  area 

vR 

Effective  rudder  inflow  velocity 

<*R 

Effective  rudder  inflow  angle 

Ur.Vr 

Components  of  rudder  effective  inflow  velocity 

€ 

constant  in  Equation  (2.7c) 

UP 

Effective  propeller  inflow  velocity 

k 

constant  in  Equation  (2.7c) 

JVp 

Thrust  coefficient 

J 

Advance  coefficient 

n 

Number  of  propeller  revolutions  per  second 

D 

Propeller  diameter 

wp 

Effective  propeller  wake  fraction 

T 

constant  in  Equation  (2.7e) 

*P 

x  coordinate  of  propeller  position 

Cpv»Cpr 

propeller  flow  rectification  coefficients 

Y 

flow  rectification  coefficients 

CSr'C5rrr»CSnv 

rudder  wake  coefficients 
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FN 


(  6A3A\ 
(a+2.25J 


(A  \ 

^R 


\P 


V*  sinCa,) 


V*    -    f». 


2  2 

+    V 
R  VR 


(2.7a) 
(2.7b) 


=    V 


\ 


1  + 


(2.7c) 


J    =    -*- 
nD 


Up    =   "[(l-Wp+T^V+^+C^V+C^}] 


aR    =    6     +    tan"1 


/v  \ 


(2.7d) 

(2.7e) 
(2.7f) 

(2.7g) 


E.      RIGHTING  MOMENT  REPRESENTATION 

For  the  initial  part  of  the  righting  arm  curve  and  for  wall  sided  ships,  Equation 
(2.8)  gives  a  good  approximation  for  the  function  GZ(<£).  [Ref.  16] 


GZ(4>)    =    GMsin(<t>)     + 


BM 


tan2((j))sin(<J>) 


(2.8) 


GM    =    the  transverse  metacentric  height 


where 


BM    =    the  transverse  metacentric  radius 
<j)  =  the  roll  angle 


Equation  (2.8)  expresses  the  fact  that  for  most  ships  of  fairly  rectangular 
midship  section,  the  GZ(<f>)  curve  exhibits  the  typical  characteristics  of  a  hardening 
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spring  restoring  moment.  The  spring  constant  is  an  increasing  function  of  the  roll 
angle  </>.  To  see  this,  we  use  the  Taylor  series  expansion  up  to  third  order  for  the 
trigonometric  functions: 

sin(4>)    =    <J>    -    -  <t>3 
o 

tan(<|>)    =    <|>    +    |  e 
Therefore,  if  terms  higher  than  third  order  are  neglected: 


-H 


W  1       \  /  1       \ 

31  (a      1a3|    _    x2[x       *x3 


tanz(<|>)sin(<|))  =  |4>  +  t<M  |4>-g4n  =  r|4>"T* 


and  substituting  into  Equation  (2.8)  yields: 


-   A3 


=  **, 


GZ(<J>)  =  gm(<|>--4>3' 


1^,*  -^Ti,  (l^Tlr  1 


+  -#M(|>3  =  GA/<|>  +  l-BM      -GM 


4>3 


Now  for  most  ships,32JA/>GA/>  which  means  that  the  leading  cubic  coefficient 
the  GZ($)  curve  is  positive  and  the  curve  concaves  up  as  would  a  hardening 


m 
spring. 


F.       STATE-SPACE  REPRESENTATION 

Combining  Equations  (2.2)  and  (2.4)  and  rearranging  the  terms,  the  equations 
of  motion  can  be  expressed  in  vector  form: 
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Z     =    F(z,6) 


(2.9) 


where         z    is  the  state  vector    [v  r  p  <|>]r. 
Equation  (2.9)  is  more  conveniently  expressed  in  matrix  notation  as  is  shown  in 
Equations  (2.10).  The  last  matrix  equation  is  derived  from  the  kinematic  roll  rate. 


{mxG-NJ     (IrNJ         0        -ty 
(m-y.)     (mxG-Yf)     -mzG     -Y± 


-mz, 


0 


0         (/,-*,)    -% 
0  0  1 


V 

fr 

r 

fs 

P 

fR 

4>. 

P. 

fr    =    Nvv+N^v3+N„VT2+(Nr-mxGu)r+Nrr/*+Nvvrv2r+N^  + 
^y^N^v^N^r^N^r^N^)^ 

fs    =    nv+lrWvV3^VTrvr2+(rr-m«)r+rfr/3+ywv2r+y+<t)+ 
Y^+Y^+Y^+Y^rtf+YW+Yc 

fR    =    ^vV+ATwv3^wv2r+(^+mzG«)r+/i:frrr3+^vr2+^<|)  + 
K^v^K^v^^K^r^K^r^Kib^K.-AGZ^) 


(2.10a) 


(2.10b) 


(2.10c) 


(2.10d) 
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III.   STEADY  STATE  REPRESENTATION 

A.      LINEAR  STEADY  STATE  SOLUTIONS 

With  the  above  third  order  Taylor  expansion,  the  equations  of  motion 
represented  by  Equations  (2.10)  are  highly  nonlinear,  the  solution  of  which  depends 
upon  the  initial  and  equilibrium  conditions.  The  steady  state  condition  requires  the 
acceleration  terms  Vj  f  a   ^  to  be  equal  to  zero.  Thus  Equation  (2.10a)  reduces  to 

Equation  (3.1)  for  steady  state: 


ft 

fs 

f, 

p\ 


-    [0]  (3-1) 


where  fY,  fs,  and  fR  are  defined  in  Equations  (2.10). 

The  numerical  solution  to  Equation  (3.1)  is  an  iterative  process  which  depends 
upon  a  fairly  accurate  initial  estimate  of  the  root  for  convergence.  This  starting  point 
can  be  determined  by  assuming  an  equilibrium  state  of  v0=r0=p0=</>0=0  and  a  very 

small  initial  rudder  deflection  ( 60£l°).  Application  of  these  assumptions  to  Equation 

(3.1)  yields  the  linear,  steady-state,  yaw-sway-roll  equations  of  motion  in  Equation 
(3.2). 
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Nv  (Nrmx^)         N+ 

V 

-N(6)-N0 

Yv     (Yr-mu)            ^ 

r 

= 

-Y(6)-YQ 

Ky  (Kr+mZGu)   (K+-AGM) 

.*. 

-K(6)-K0 

(3.2) 


Equation  (3.2)  is  in  the  familiar  matrix  form  Az  =  B,  where  z  is  the  state  vector 
with  p=0. 

Various  computer  subroutines  are  commercially  available  which  solve  this 
simple  system  of  equations  for  the  state  vector  values,  given  the  other  parameters. 
The  values  for  the  state  vector  obtained  in  this  manner  are  then  used  as  the  initial 
estimate  in  solving  the  nonlinear  system  of  equations. 

B.       NONLINEAR  STEADY  STATE  SOLUTION 

Once  the  initial  estimate  for  equilibrium  is  obtained,  the  system  of  nonlinear 
equations  in  Equation  (3.1)  can  be  solved  numerically  for  the  state  vector  nominal 
points  for  each  rudder  angle,  8.  Here  again,  commercial  computer  subroutines  are 
available.  The  Levenberg-Marquardt  algorithm  with  analytic  Jacobian  was  chosen  as 
the  solution  method  [Ref.  17].  Since  p=0  in  Equation  (3.1),  the  system  of 
nonlinear  equations  is  reduced  to  three  and  the  Jacobian  was  determined  by 
Equations  (3.3). 
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M    = 


dfy  dfy  dfy 

dv  dr  dfy 

%  %  % 

dv  dr  dfy 

%  %  % 

dv  dr  dfy 


(3.3a) 


* 


-^  =  Nv+3N„v2+N^2+2N„vr+2N^+N^+2!!&l 


dv 


dv 


^=2N^vrHNrmx^y3N^2+N^v2+2N^r^N^2JN(b} 


dr 


dfy 


S  =  ^+^v2+2^v*+^r2+2^r<D  +  ^> 


J  =2ywvr+(rr-m«)+3y^2+rwv2+2rrr,r(|>+rnH4>2+-^ 
§  =  V^v2+2y^v*+y^r2+2y^r<|>+^ 

^  =  ^v+3^vvvV2^^2+2^vr+2^v*+^^2+-?|P 

^=2^vr+(^r+m2^)+3/r^2+^v2+2^r*+/r^2+^ 

§  =^+^v2+2^^+^2+2^*+^-M^ 


dfy 


dfy  dfy 


(3.3b) 
(3.3c) 

(3.3d) 
(3.3e) 

(3-3f) 

(3.3g) 

(3.3h) 
(3-3i) 

(3-3j) 
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tf(6)    =    F,(v^) 

w/^re  1TO    =    F2(v,r) 

tf(6)    =    F3(v/) 


and 


a(AGZ(4>))  _  A 


dN(6) 

a<t> 


=  o 


=  o 


=    0 


GA/cos(4>)+BA/tan2(<t>) 


{cos($)    2 


cos(<J>) 


(3.3k) 


(3.31) 


Since  Fx,  F2,  and  F3  are  complicated  functions  of  v  and  r  (Equations  (2.6)  and 
(2.7)),  the  following  approximation  method  was  used  for  evaluation: 


-1      I  equilibrium  =x„ 


FjxJ  -  F{0.99x<) 


(3.4) 


Each  term  in  the  Jacobian  is  determined  by  evaluation  at  the  previous  nominal 
point  with  the  initial  nominal  point  determined  from  linear  solutions  as  described  in 
Section  A  above. 

The  result  of  the  numerical  solution  to  the  nonlinear  equations  (Equations 
(3.1))  is  the  steady  state  values  for  v,  r,  and  <£  for  each  rudder  angle  under  steady 
turning  conditions  and  for  a  given  GM  and  propeller  speed. 

C.      TYPICAL  RESULTS 

Utilizing  the  solution  method  depicted  in  Sections  A  and  B  above,  PROGRAM 
COUPLED  listed  in  Appendix  A,  was  developed  to  predict  the  steady  state  roll 
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angle,  sway  velocity  and  turning  rate  of  a  high  speed  container  ship  as  a  function  of 
rudder  angle. 

The  design  characteristics  of  a  SR-108  container  ship  are  presented  in  Table 
IV.  The  hydrodynamic  derivatives  and  coefficients  used  in  the  computer  simulation 
are  listed  in  Table  V. 

For  a  base-line  model,  the  hydrodynamic  derivatives  and  coefficients  were  held 
constant  in  addition  to  the  Froude  number  (hence  the  propeller  revolutions),  the 
transverse  metacentric  height,  and  speed  loss  ratio.  Therefore,  for  the  figures 
presented  in  this  section: 

Froude  number  Fn  =  0.3 

Metacentric  height  GM  =  0.3  m 

Speed  loss  ratio  a  =  0.6 

Figures  3.1  through  3.3  show  the  variation  of  sway  velocity,  yaw  rate,  and  roll 
angle  as  a  function  of  rudder  angle  during  a  steady  turn.  It  can  be  seen  that  all 
steady  state  variables  v,  r  and  <f>  are  highly  nonlinear  functions  in  8  which  is 
attributed  to  both  the  effect  of  the  nonlinear  terms  in  our  equations  of  motion  and 
the  fact  that  the  ship  speed  is  reduced  during  the  turn.  Had  a  constant  speed,  linear 
model  been  utilized,  v,  r  and  </>  would  appear  as  straight  lines  versus  S. 
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TABLE  IV.   PRINCIPAL  DIMENSIONS  OF  SR-108  CONTAINER  SHIP  [Ref.  9] 


Items 

Ship 

Model 

Hull  Length  B.  P. 

L 

(m) 

175.00 

3.000000 

Breadth 

B 

(m) 

25.40 

0.435000 

Draft  Fore 

TF 

(m) 

8.00 

0.137100 

Aft 

TA 

(m) 

9.00 

0.154300 

Mean 

T 

(m3) 

8.50 

0.145700 

Displacement  Volume 

(m3) 

21,222 

0.106860 

Height  from  keel  to 

transverse  metacenter 

KM 

(m) 

10.3900 

0.17810 

Height  from  keel  to 

center  of  buoyancy 

BM 

(m) 

4.6154 

0.07912 

Block  coefficient 

Cb 

0.55900 

Prismatic  coefficient 

Cb 

0.58000 

Waterplane  area 

0.68600 

coefficient 

Cvv 

0.518  L 

Midship  section 

Cm 

0.240  L 

coefficient 

LCB  from  forward 

Perpendicular 

Bilge  keel 

Length 

Ld 

(m) 

43.75 

0.7500 

Depth 

(cm) 

45.00 

0.7714 

Rudder  Area 

AR 

(m2) 

33.0376 

0.009709 

Height 

H 

(m) 

7.7583 

0.133000 

Aspect  Ratio 

A 

1.821900 

Area  Ratio 

AR/Ld 

(m) 

1/45.0 

Propeller 

Diameter 

D 

(m) 

6.533 

0.112 

Pitch  Ratio 

P 

1.009 

Expanded  Area  Ratio 

0.670 

Boss  Ratio 

0.180 

Number  of  blades 

5 

24 


TABLE  V.    HYDRODYNAMIC  DERIVATIVES  AND  COEFFICIENTS  [Ref.  9] 


A.)   Hull  Only 


Yv 

-0.0116000 

Nv 

-0.0038545 

K 

0.00030260 

Yr 

0.0024200 

Nr 

-0.0022200 

K, 

-0.00006300 

Y, 

-0.0000630 

N, 

-0.0001424 

K* 

-0.00002100 

Y 

-0.1090000 

N™ 

0.0014920 

Kwv 

0.00284300 

Ym 

-0.0405000 

Nm 

0.0015600 

K. 

0.00105650 

Ym 

0.0017700 

Nm 

-0.0022900 

K. 

-0.00004620 

Y 

0.0214000 

Nw, 

-0.0424000 

Kwr 

-0.00055800 

Y 

0.0460500 

NW(£ 

-0.0190580 

^w^ 

-0.00120120 

Yv<^ 

0.0030400 

Nv<w, 

-0.0053766 

Kv^ 

-0.00007930 

*rr<£ 

0.0093250 

N* 

-0.0038592 

*Nr* 

-0.00024300 

^r^ 

-0.0013560 

Nr^ 

0.0024195 

Kr** 

0.00003569 

B.)   Propeller  and  Rudder 


aH 

0.237 

T 

1.090 

(l-wp) 

0.816 

xH 

-0.480 

6 

0.921 

Xr 

-0.500 

CRX 

0.710 

k 

0.631 

*P 

-0.526 

Zr 

0.033 

c6r 

-0.156 

.K.J 

0.527-0.455J 

CPV 

0.000 

C5rrr 

-0.275 

y 

0.088  v>0 
0.193  v<0 

Cpr 

0.000 

CSrrv 

1.960 

Np 

79.10  Fn=0.2 
118.64  Fn=0.3 
158.19  Fn=0.4 
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Figure  3.1:   Sway  Velocity  Under  Steady  Turning 


Figure  3.2:   Turning  Rate  Under  Steady  Turning 
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Figure  3.3:    Roll  Angle  Under  Steady  Turning 
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IV.   STABILITY  ANALYSIS 


A.      COUPLED  STEERING  AND  ROLL  EQUATIONS  OF  MOTION 

Once  the  nominal  points  are  determined  for  each  rudder  angle,  Equation 
(2.10a)  can  be  reduced  to  state-space  representation  (Equations  (3.4))  for  stability 
analysis. 


where 


z    = 


B    = 


Bz    =    Az 

v 
f 
P 


z    = 


v 

r 

P 

4>J 


(mxG-NJ     (7Z-Ay  0        -ty 

(m-YJ     (mxG-Yf)     -mzG     -Y± 


mz, 


0 


A    = 


0 
0 

dfy  dfy 

dv      dr 

%     Qfs 
dv      dr 

%     % 
dv      dr 

dp      dp 
dv      dr 


o 


K, 


dfy 

dfy 

~dp 

3<t> 

% 

Vs 

dp 

&t> 

% 

% 

dp 

<3<f> 

dp 

Jp 

dp 

a4>. 

(4.1a) 


(4.1b) 


(4.1c) 
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Matrix  B  can  be  thought  of  as  the  generalized  mass  matrix  for  the  problem 
while  matrix  A  is  the  Jacobian  matrix  of  f(z,S)  evaluated  at  each  nominal  point.  For 
stability  analysis  purposes,  the  variables  v,  r,  p,  and  <f>  are  understood  to  represent 
small  deviations  of  the  actual  variables  from  their  respective  steady  state  values  as 
computed  in  the  previous  chapter.  The  rationale  behind  Equations  (4.1)  is  the 
introduction  of  small  local  perturbations  in  v,  r,  p,  and  (f>  superimposed  on  their 
steady  turning  values.  Lyapunov's  linearization  theorem  establishes  that  stability 
properties  of  a  nominal  point  for  a  nonlinear  system  can,  in  general,  be  deduced 
from  the  stability  properties  of  the  corresponding  linearized  system. 

The  terms  in  the  matrix  A  are  as  defined  in  Equations  (3.3)  with  the  following 
additions: 


dfy 

Vs 

Vr 

r\ 

=     0 

dp 

dp 

dp 

dp     . 

.     dp     _ 

dp 

=    0 

dv 

dr 

d<t> 

dp     _ 

1 

dp 

Stability  can  be  determined  from  the  solution  to  the  eigenvalue  problem: 

XBz    =    Az 

(4.2) 

L4  -  XB\     =0 
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If  all  eigenvalues,  A,  have  negative  real  parts,  the  nominal  point  is  asymptotically 
stable.  If  at  least  one  eigenvalue  is  positive,  the  nominal  point  is  unstable.  The 
degree  of  stability  of  the  nominal  point  is  therefore  related  to  the  real  part  of  the 
eigenvalues.  The  more  negative  they  are,  the  faster  the  exponential  convergence  of 
solutions  in  time  to  the  nominal  point.  The  imaginary  part  of  the  eigenvalues 
characterizes  the  frequency  content  of  the  oscillatory  behavior  of  the  system 
response. 

A  subroutine  was  included  in  the  computer  model  to  solve  this  eigenvalue 
problem  as  a  function  of  rudder  angle. 


B.      DECOUPLED  STEERING  EQUATIONS 

In  order  to  distinguish  the  characteristic  contributions  of  the  sway  and  yaw 
motions  from  those  of  roll,  the  steering  equations  (Equations  (4.3))  were  decoupled 
from  the  maneuvering  equations  (Equations  (4.1)),  resulting  in  Equations  (4.3). 


(mxG-NJ      (Iz-Nr) 
(m-YJ      imxG-Yr) 


Jys 
fss 


(4.3a) 


frs  =  Ys+Y^+Y^+iYr-muy+Y^+Y^r+Yib)  (4-3b) 

frs  =  N^N^vKN^r2^Nr-mxGu)r^rr/^Nwrv1r^N{b)  (43c) 

This  decoupling  is  based  on  the  assumption  that  all  cross  coupling  coefficients 
between  roll  and  sway/yaw  are  zero,  which  is  the  usual  approximation  made  for 
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surface  ships  where  roll  is  studied  independently  from  sway/yaw.  One  of  the 
motivations  for  this  study  is  to  evaluate  the  degree  of  accuracy  of  such  an 
approximation. 

The  additional  subscript,  S,  distinguishes  the  decoupled  steering  forces  and 
moments  from  those  of  the  fully  coupled  system.  Following  the  same  format  as 
before,  the  state-space  formulation  for  stability  analysis  is  given  in  Equations  (4.4). 


where 


B< 


Bsi&  = 

A-s^z 

V 

V 

Ss  ■ 

r 

t 

5s  = 

r 

(mxG-N,)  (Iz-Nr) 

(m-YJ  (mxG-Y,) 

dv  dr 

&ss  &ss 

dv  dr 


**     =    N+3Nv2+Nr2+2Nvr+dNib) 


dv 


v  vw  vrr  wr 


dv 


dr 

Vss 

dv 

&ss 
dr 


=  2Nvr+(N-mxru)+3Nr2+Nv2+-dN^ 

vrr         v     r  O    '  rrr  wr 


=    Y   +  3Y    v2  +Yr2+2Y   vr 

v  vw  vrr  vrr 


dr 

dY(6) 

dv 


=  2rvrrvr+(yr-m«)+3yrrrr2+ywv2+-^ 

or 


(4.4a) 


(4.4b) 


(4.4c) 


(4.4d) 
(4.4e) 
(4.4f) 
(4-4g) 
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The  approximation  method  of  Equation  (3.4)  was  used  to  determine  rudder 
force  and  moment  contributions  in  the  Jacobian  matrix  A. 

The  stability  attributable  to  yaw  and  sway  as  a  function  of  rudder  angle  is  thus 
determined  by  the  solution  of  the  decoupled  steering  eigenvalue  problem  in 
Equations  (4.5). 


*s~s£s         AsZf 


or 

\AS-XSBS\     =    0 


(4.5) 


C.      DECOUPLED  ROLL  EQUATIONS 

By  decoupling  the  steering  equations,  the  roll  equations  become  uncoupled 
producing  the  following  relationship: 


0  1 


JRR 
P 


(4.6a) 


fm    =    KA  +  K(b)  -  AGZ(d>) 


(4.6b) 


The  additional  subscript,  R,  refers  to  the  decoupled  roll  forces  and  moments 
and  GZ((f>)  is  as  given  in  Equation  (2.8).  Equations  (4.7)  are  the  state-space  stability 
analysis  formulations. 
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Br*e    ~    ar^r 


where 


ht 


>     £i 


(4.7a) 


B„    = 


0 


-*. 


(4.7b) 


%R 

#** 

a? 

3p 

^> 

3p 

L41 

34, 

(4.7c) 


« 


W? 


dp 


=    0 


%r    =    K     a(AGZ(4))) 

a<j>         ♦       64> 


(4.7d) 


3p    =   1 

c3p 


Sl  =  o 

34> 


The  expression  — ^!2Z  is  given  in  Equation  (3.3). 

d4> 


(4.7e) 


Therefore,  the  stability  attributable  to  roll  as  a  function  of  rudder  angle  is 
attained  by  solution  of  the  decoupled  roll  eigenvalue  problem  of  Equations  (4.8). 


*-rBrZr    "    ^rIr 


or 


(4.8) 
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D.      TYPICAL  RESULTS 

Continuing  with  the  base-line  model  described  in  Chapter  III  Section  C,  the 
roots  of  Equations  (4.2),  (4.5)  and  (4.8)  were  obtained  using  PROGRAM 
COUPLED  (Appendix  A).  Figure  4.1  is  a  plot  of  the  real  roots  of  the  fully  coupled 
maneuvering  equations  as  a  function  of  rudder  angle.  As  can  be  seen,  there  are  three 
distinct  components  which  were  classified  as  due  to  roll  or  due  to  steering  (sway/yaw) 
by  decoupling  the  steering  and  roll  equations  as  discussed  in  Sections  B  and  C  above. 
Figures  4.2  through  4.4  show  the  three  components  of  the  fully  coupled  system 
matched  with  their  respective  decoupled  solutions  of  Equations  (4.5)  for  steering  and 
Equations  (4.8)  for  roll. 

In  Figure  4.2,  the  decoupled  roll  real  roots  are  constant  for  each  rudder  angle, 
whereas  the  coupled  roots  vary  with  an  increasing  stability  to  about  five  degrees  of 
rudder.  The  upper  steering  component  shown  in  Figure  4.3  indicates  very  little 
deviation  in  the  coupled  and  decoupled  roots,  however,  the  lower  steering 
component  in  Figure  4.4  shows  a  deviation  in  the  two  roots  as  the  rudder  angle  is 
increased.  The  steering  eigenvalues  become  increasingly  more  negative  as  the  rudder 
angle  is  increased.  This  indicates  that  the  steering  system  is  dynamically  more  stable 
for  a  non-zero  rudder  angle,  and  as  a  result,  a  steering  system  that  is  designed  to  be 
stable  for  straight  line  motions  will  be  even  more  stable  for  motions  along  curved 
reference  paths. 

The  imaginary  roots  for  the  decoupled  equations  are  zero,  indicating  that  all 
imaginary  roots  obtained  in  the  solution  of  the  coupled  equations  are  due  to  roll 
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only.  Figure  4.5  shows  that  for  rudder  angles  less  than  15°,  the  coupled  model 
produces  larger  values  than  the  decoupled  model  and  for  rudder  angles  greater  than 
15°,  the  reverse  is  true.  Comparing  Figures  4.2  and  4.5,  we  can  see  that  the  imaginary 
part  of  roll  response,  or  the  roll  frequency  of  oscillation,  can  be  predicted  fairly 
accurately  from  the  decoupled  roll  equations.  This  is  not  the  case  for  the  real  part 
where  decoupling  the  equations  results  in  a  severe  underestimation. 

Figure  4.6  is  a  root  locus  plot  of  the  roll  eigenvalues  for  both  coupled  and 
uncoupled  models.  As  can  be  seen,  the  fully  coupled  model  is  more  stable  than  the 
decoupled  model. 

The  roll  damping  caused  by  roll  is  shown  in  Figure  4.7.  The  decoupled  roll 
analysis  shows  less  damping  for  each  rudder  angle  than  does  the  fully  coupled  model. 
Furthermore,  the  damping  increases  significantly  for  large  rudder  angles.  It  can  be 
seen  that  the  actual  roll  damping  is  significantly  larger  than  what  is  predicted  from 
the  decoupled  model,  which  may  underestimate  the  actual  damping  ratio  by  a  factor 
of  two  to  three. 
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Figure  4.1:   Coupled  Model  Real  Eigenvalues 
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Coupled  Model  and  Decoupled  Model  Upper  Steering  Eigenvalues 
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V.   SENSITIVITY  OF  MATHEMATICAL  MODEL 
TO  VARIATION  IN  DESIGN  PARAMETERS 


A.   BACKGROUND 

Since  actual  ships  rarely  operate  under  the  precise  conditions  assumed  in  the 
development  of  this  mathematical  model,  changes  in  certain  variables  to  measure  the 
response  sensitivity  is  prudent,  as  for  any  design  study.  The  following  sensitivity 
analysis  has  this  aim  in  mind  as  well  as  to  identify  any  inexplicable  response  as  an 
area  for  further  study.  It  must  be  remembered  that  certain  assumptions  were  made 
at  the  outset  which  might  need  refinement  as  the  design  spiral  continues.  For 
instance,  the  roll  damping  term,  K^,  is  assumed  to  be  constant  throughout  the 
simulation,  although  experimental  evidence  [Ref.  18]  suggests  that  roll  damping 
coefficients  are  forward  speed  and  motion  amplitude  dependent.  By  isolating  a 
parameter,  problem  areas  are  enhanced  so  that  a  particular  parameter  can  be  the 
focus  for  analysis  of  the  validity  of  the  assumptions  made  or  the  modelling  method 
used. 

If,  on  the  other  hand,  one  or  more  of  the  dependent  variables  do  not  show 
appreciable  sensitivity,  further  simplifying  assumptions  could  possibly  be  made  when 
analyzing  the  more  difficult  problems  of  maneuvering  in  regular  and  irregular  seas. 

In  the  following  sections,  sensitivity  to  speed  loss  during  a  turn,  Froude  number 
of  approach  spped,  metacentric  height,  operating  turn,  and  righting  moment  is 
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evaluated.  For  brevity,  a  limited  number  of  plots  are  presented  in  each  case. 
Appendix  B  contains  additional  plots  for  metacentric  height,  Froude  number,  and 
speed  loss. 

B.       SPEED  LOSS  DURING  A  TURN  AS  THE  DESIGN  PARAMETER 

The  sensitivity  of  the  model  to  changes  in  speed  loss  during  the  turn  was 
investigated  by  holding  transverse  metacentric  height  (GM)  and  Froude  number 
constant: 


GM    =    0.3  m 
Fn    =    0.3 

In  determining  the  speed  loss  under  steady  turning  conditions,  it  was  assumed 
that  the  water  depth  was  sufficient,  i.e.  >  110  meters,  to  neglect  this  effect  [Refs.  2, 
3].  However,  it  has  been  determined  by  several  sources:  Davidson  (1944),  Shiba 
(1960),  Strom-Tejsen  (1965),  Eda  and  Crane  (1965),  among  others  [Ref.  2],  that 
speed  loss  in  steady  turning  is  a  function  of  the  hull  configuration  (including  rudder) 
and  the  turning  diameter.  Various  other  sources,  [Refs.  3,  9,  10]  consider  rudder 
angle  as  an  additional  parameter  in  determining  speed  loss.  In  an  effort  to  combine 
these  effects,  Figures  25  and  178  and  Tables  29(a)  and  29(b)  from  Principles  of  Naval 
Architecture  [Ref.  2],  Figure  11  of  Eda  and  Falls  [Ref.  3],  and  Figure  10  of  Son  and 
Nomoto  [Ref.  9]  were  used  to  develop  the  following  simple  relationship  to  represent 
speed  loss: 
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—     =    1  -  a  *  6r 
V 

V    =    steady  turning  speed 

V0  =    approach  speed 

where 

a     =    speed  loss  ratio 

6r  =    rudder  angle  in  radians 

Variation  of  the  speed  loss  parameter  ranges  from  no  speed  loss  (a=0)  to  maximum 

speed  loss  (a=0.73)  during  a  turn  with  fixed  controls.  Results  of  these  effects  on 

sway  velocity,  yaw  rate  and  roll  angle  may  be  viewed  in  Figures  5.1  through  5.3.  The 

influence  of  speed  loss  on  r  is  negligible,  indicating  that  in  all  cases  considered,  the 

ship  steady  turning  rate  performance  is  not  affected.  On  the  contrary,  v  and  4>  show 

a  relatively  significant  modification  due  to  speed  loss.  The  sign  reversal  of  the  roll 

angle  plot  (Figure  5.3)  for  large  rudder  angles  is  rather  uncommon,  but  it  is, 

nevertheless,  possible  as  a  result  of  the  nonlinearities  that  are  present  in  the 

equations  of  motion  for  this  model. 

The  characteristic  roots  as  a  function  of  rudder  angle  were  determined  to 

establish  stability  trends.  The  real  roots  indicate  three  distinct  components  as 

previously  discussed  in  Chapter  IV,  Section  C.  Both  the  upper  and  lower  steering 

components  demonstrate  a  significant  increase  in  stability  with  rudder  angle  but  little 

variance  in  stability  due  to  speed  loss.  The  roll  component,  however,  indicates 

greater  stability  as  speed  loss  decreases.  We  discovered  earlier,  that  the  non-zero 

values  in  the  root  locus  were  due  solely  to  the  roll  contribution  owing  to  the 

imaginary  parts  of  the  steering  elements.  The  damping  ratio  of  the  roll  contribution 
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Figure  5.1: 


Effect  of  Speed  Loss  on  Sway  Velocity 
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Figure  5.3:  Effect  of  Speed  Loss  on  Roll  Angle 


imaginary  parts  of  the  steering  elements.  The  damping  ratio  of  the  roll  contribution 
is  related  to  both  real  and  imaginary  parts  of  the  eigenvalues.  Figures  5.4  and  5.5, 
for  the  fully  coupled  and  uncoupled  models  respectively,  show  this  as  a  function  of 
rudder  angle.  Both  models  exhibit  qualitatively  the  same  trend,  with  the  uncoupled 
model  showing  a  significant  underestimation  of  the  true  damping  ratio  value. 

C.      FROUDE  NUMBER  AS  THE  DESIGN  PARAMETER 

In  testing  the  model  sensitivity  to  Froude  number  changes,  the  transverse 
metacentric  height  and  speed  loss  ratio  were  kept  constant: 

GM    =    0.3  m 
a     =    0.6 
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Figure  5.5: 


Effect  of  Speed  Loss  on  Decoupled  Model  Damping  Ratio 
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Froude  numbers  used  to  calculate  the  forward  velocity  prior  to  initiating  the  turn 
aswell  as  propeller  revolutions  were  assumed  to  remain  constant  during  the  turn.  For 
each  Froude  number  and  rudder  angle,  the  sway  velocity,  yaw  rate  and  roll  angle 
were  determined  using  the  model  program.  Figures  5.6  through  5.8  at  the  end  of  this 
section  depict  graphic  results.  Note  that  v  and  r  are  not  readily  affected  by  change 
in  Froude  number,  whereas  (f>  is  greatly  changed,  increasing  dramatically  with  Froude 
number. 

Using  these  values,  the  stability  of  the  coupled  response  was  determined.  The 
real  roots  of  the  coupled  maneuvering  equations  for  the  three  Froude  numbers  are 
indistinguishable  when  plotted  as  in  Figure  4.1,  as  are  the  roots  for  the  lower 
steering  components  (Figure  4.2).  The  upper  steering  component  exhibits  nonlinear 
behavior  through  five  degrees  of  rudder  for  all  Froude  numbers  with  increasing 
stability  as  both  Froude  number  and  rudder  angle  increase.  The  roll  component  was 
nonlinear  throughout  the  range  of  rudder  angles  tried.  Furthermore,  as  the  Froude 
number  increased,  so  did  the  range  of  rudder  angles  for  which  stability  increased. 

The  coupled  model  roll  damping  ratio  remains  relatively  insensitive  to  Froude 
number  (Figure  5.9),  while  the  uncoupled  model  suggests  a  wide  variation  (Figure 
5.10).  This  demonstrates  the  incorrect  conclusions  which  can  be  reached  when  using 
linear/uncoupled  models  that  neglect  the  actual  coupling  of  the  lateral  sway/yaw 
motions  back  into  roll. 


46 


0.25 
0.2 

'^^ 

~~| 

\ 

\ 

Fn 

0.15 

V  ■ 

0.2 

0.3    

" 

fc     o.i 

3  °05 

>        0 

\. 

0.4    -.- 

1 

|  -0.05 
*     -01 

I 

-0.15 

- 

\ 

-0.2 
-0.25, 

- 

v^_ 

0 

10 

-20 

-10                   0                    10 

20                    3 

RUDDER  ANGLE   (deg) 

Figure  5.6:  Effect  of  Fn  on  Sway  Velocity 


1.5 


g  °'5 

£  o 
< 

u  -0.5 


-1 


-1.5 


30 


— 1                                     1                                     1 

r  — 

- 

Pn 

0.2 

0.3    — 
0.4    -.- 

- 

Jr 

- 

-20 


-10  0  10 

RUDDER  ANGLE   (deg) 


20 


30 


Figure  5.7: 


Effect  of  Fn  on  Yaw  Rate 


47 


20 

10 

0 

-      .-''' 

\                            i 1 — 

I 
1 
i 

J 
| 

Pn 

0.2 

0.3    — ■ 
0.4    -.- 

00 

-^ 

■a 
UJ 

-J 
O 

d 

0 

a: 

-10 

V 

.--" 

/_ 

IE 

-20 

in 

;\V^_ 

..--—-' 

V 

-30 


-20 


-10  0  10 

RUDDER  ANGLE   (deg) 


20 


30 


Figure  5.8:  Effect  of  Fn  on  Roll  Angle 
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Figure  5.9:  Effect  of  Fn  on  Coupled  Model  Damping  Ratio 
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Figure  5.10:  Effect  of  Fn  on  Decoupled  Model  Damping  Ratio 

D.      TRANSVERSE  METACENTRIC  HEIGHT  AS  THE  DESIGN  PARAMETER 

For  transverse  metacentric  height  as  the  design  parameter  in  the  model,  Froude 
number  and  speed  loss  ratio  were  held  constant: 

Fn    =    0.3 
a     =    0.6 

Here  again,  the  sway  velocity,  yaw  rate  and  roll  angle  for  each  GM  and  rudder 
angle  were  computed  using  the  model  program.  The  results  are  represented  in 
Figures  5.11  through  5.13.  As  indicated  for  the  case  of  Froude  number  being  the 
design  parameter,  v  and  r  are  affected  very  little  by  even  an  order  of  magnitude 
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change  in  GM.  On  the  other  hand,  the  increase  in  GM  not  only  causes  a  large 
decrease  in  the  magnitude  of  the  roll  angle,  but  in  the  range  as  well. 

In  determining  the  influence  of  GM  on  the  stability  characteristics,  it  was  found 
that  the  real  roots  plotted  as  a  function  of  rudder  angle  showed  no  distinction 
between  GM's.  When  the  components  were  separated  out,  only  the  upper  steering 
and  roll  constituents  exhibited  significant  variation.  The  upper  steering  response 
reflected  a  somewhat  incongruous  reaction,  that  is,  as  the  GM  increases,  the  stability 
due  to  steering  decreases  slightly.  The  roll  component  demonstrated  a  marked 
increase  in  stability  with  increase  in  GM.  This  result  is  evident  in  the  root  locus 
diagram  of  Figure  5.14. 

Conversely,  the  damping  ratio  did  not  change  significantly,  as  shown  in  Figure 
5.15.  Computation  of  the  damping  ratio  from  the  decoupled  model  showed  a  much 
wider  deviation  among  different  GM's. 

E.      TRIM  AS  THE  DESIGN  PARAMETER 

For  the  investigation  of  a  vessel  in  the  trimmed  condition,  the  transverse 
metacentric  height,  Froude  number  and  speed  loss  ratio  were  held  constant: 


GM    =    0.3  m 
Fn    =    0.3 
a     =    0.6 


50 


0.3 


0.2- 


H 

01 

U 

3 

u 

> 

0 

>- 

< 

£ 

CO 

-0.1 

-0.2- 


•°-?30 


-20 


-10  0  10 

RUDDER  ANGLE   (deg) 


ca*  (■) 

0.3    

1.0 

3.0    -.- 


20 


30 


Figure  5.11:  Effect  of  GM  on  Sway  Velocity 
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Figure  5.12:  Effect  of  GM  on  Yaw  Rate 
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Figure  5.13:         Effect  of  GM  on  Roll  Angle 
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Figure  5.15:  Effect  of  GM  on  Coupled  Model  Damping  Ratio 
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Figure  5.16:         Effect  of  GM  on  Decoupled  Model  Damping  Ratio 
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The  base-line  model  chosen  for  this  study,  and  therefore  the  hydrodynamic 
derivatives  and  coefficients  derived  from  model  tests,  were  for  a  design  trim  of  0.5 
meters  by  the  stern  and  a  mean  draft  of  8.5  meters  for  the  actual  ship  (refer  to  Table 
IV  for  model  parameters).  It  has  been  shown  [Ref.  11]  through  linear  theory  that 
these  derivatives  and  coefficients  change  for  the  trimmed  condition  by  the  following 
semi-empirical  relations:  [Ref.  15] 


w 

-  rv«»[i  ♦  32;J 

Nv(x) 

■  tfv(0) 

0.27 
1  -  

Id 

v    m 
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=    Yf(0) 

i  +  °-8* 

N,W 

=    Nr(0) 

i  ♦  °-3* 

K  = 

Yv(0) 

where  +   for  fr/m  ^V  stern 

T       = 

-   for  trim  by  bow 
dm    =    design  mean  draft 

YV(Q)    =    Y  at  design  condition 


The  application  of  these  relationships  are  somewhat  questionable,  but  in  the 
absence  of  model  tests,  they  will  suffice  to  indicate  trends  for  further  study. 

A  further  modification  in  the  model  was  made  to  the  zR  parameter  (Table  III) 
or  the  z  coordinate  of  the  point  on  which  the  rudder  force  acts: 
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z*W    =    zR(0)  +  A 
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Figures  5.17  through  5.19  depict  the  effects  of  —  =  ±0.2  which  equates  to  a 
significant  change  in  operating  condition.  The  sensitivity  of  v  and  r  to  the  previous 

parameters  was  practically  negligible  with  the  exception  of  the  speed  loss  ratio.  This 

is  not  the  case  for  trim  where  turning  rate  is  improved  for  trim  by  the  bow  and 

decreased  for  trim  by  the  stern.  It  should  be  noted  that  this  is  contrary  to  linear 

theory  results.  For  sway  velocity  at  small  angles  of  rudder  the  variation  is  small,  but 

for  large  rudder  angles  the  trim  by  stern  curve  deviates  considerably  from  the  other 

two  conditions.  The  difference  is  more  pronounced  in  the  roll  angle  curves  (Figure 

5.19)  where  trim  by  the  bow  generates  higher  angles  during  the  turn.  On  the  other 

hand,  trim  by  the  stern  results  in  smaller  roll  angles  with  a  roll  angle  sign  reversal 

beyond  a  certain  rudder  angle.  This  is  due  to  the  fact  that  for  large  rudder  angles, 

the  hydrodynamic  moment  exerted  by  the  rudder  force  dominates  over  the  fluid 

forces  on  the  hull. 

The  root  locus  diagram  (Figure  5.20)  for  the  roll  stability  analysis  is  revealing 

in  that  very  different  trends  can  be  observed.  Trim  by  the  stern  is  considerably  more 

stable  than  by  the  bow  for  the  entire  range  of  heel  angle,  although  at  larger  rudder 

angles,  the  stability  is  less  than  at  the  design  condition.  This  is  possibly  qualitatively 

explained  by  noticing  in  Figures  5.21  and  5.22,  that  the  damping  ratio  for  trim  by  the 

stern  increases  sharply  precisely  at  the  inflection  point  of  the  roll  angle  as  shown  in 
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Figure  5.19.  Here  again,  it  must  be  pointed  out  that  the  damping  ratio's  variance 
between  the  coupled  model  and  the  decoupled  model  is  at  least  a  factor  of  two.  It 
can  be  seen  from  both  coupled  and  decoupled  models  that  at  the  inflection  point  of 
the  roll  angle,  the  damping  ratio  reaches  its  peak  value.  This  maximum  value  of  the 
damping  ratio  exists  for  the  other  trim  conditions  as  well,  but  its  actual  location  is 
shifted  outside  the  range  of  rudder  angle  variations  considered  here. 

F.       RIGHTING  ARM  CURVE  AS  THE  DESIGN  PARAMETER 

The  sensitivity  of  the  model  to  changes  in  the  righting  arm  as  a  function  of  the 
roll  angle  was  tested  by  holding  the  transverse  metacentric  height,  Froude  number 
and  speed  loss  constant: 


GM    =    1.0  m 
Fn    =    0.3 

a     =    0.6 


In  general,  the  righting  arm  curve,  GZ(<£),  is  taken  to  be  a  very  simplified 
version  of  Equation  (2.8),  namely  GZ(</>)  =  GM*<£.  This  linear  righting  arm  gives 
good  approximations  for  small  angles  of  heel  and  is  derived  much  the  same  as  was 
done  for  Equation  (2.8)  in  Chapter  II,  Section  E.  For  the  purpose  of  this 
investigation,  the  simplified  GZ(<f>)  was  left  in  its  nonlinear  form, 
GZ(4>)  =  GM* sin ((f)).  Figure  5.23  illustrates  the  two  GZ(#)  curves  used.  Note  that 
GZ1  (representing  Equation  (2.8))  develops  a  much  larger  righting  arm  at  increased 
angles  of  heel  than  does  GZ2  (representing  GM*sin($)). 
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Figure  5.17:         Effect  of  Trim  on  Sway  Velocity 
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Figure  5.21:  Effect  of  Trim  on  Coupled  Model  Damping  Ratio 
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Figure  5.22:  Effect  of  Trim  on  Decoupled  Model  Damping  Ratio 


59 


The  effects  of  these  different  righting  arm  curves  on  the  model  for  v,  r,  and  (f>, 
are  shown  in  Figures  5.24  through  5.26.  Sway  velocity  and  yaw  rate  appear  to  be 
unaffected  by  this  parameter  change  in  contrast  to  roll  angle,  which  shows  relatively 
high  sensitivity,  as  expected. 

Figures  5.27  (coupled  model)  and  5.28  (decoupled  model)  indicate  the  same 
trend  in  damping  ratio  sensitivity  as  the  previous  parameters,  specifically,  the 
decoupled  model  underestimates  c  by  a  factor  of  two. 

G.      SUMMARY  OF  RESULTS 

A  qualitative  summary  of  results  is  presented  in  Table  VI  below.  All 
comparisons  are  to  the  base-line  configuration  of:1 


GM    =    0.3  m 
Fn    =    0.3 
a     =    0.6 


1  For  GZ  sensitivity  GM  =  1.0  m. 
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Figure  5.23:  The  Two  Righting  Arm  Curves 
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Figure  5.26:  Effect  of  GZ(tf>)  on  Roll  Angle 
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Figure  5.27:  Effect  of  GZ($)  on  Coupled  Model  Damping  Ratio 
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Figure  5.28:  Effect  of  GZ(<£)  on  Decoupled  Model  Damping  Ratio 
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Notation  in  the  table  refers  to  the  largest  noted  deviations  from  the  model 
according  to  the  following  subjective  scale: 


Sensitivity 

Symbol 

Standard 

Negligible 

N 

<  1% 

Low 

L 

1%  -  5% 

Moderate 

M 

5%  -  20% 

High 

H 

>  20% 

Mixed 

X 

=>  large  change 
over  the  range 

Additionally,  A,  B  and  C  refer  to  mdder  angle  range  for  the  noted  deviations  such 
that: 

0'  <  A  z  10° 
10°  <  B  z  20° 
20°    <    C    <l    30" 
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TABLE  VI.   QUALITATIVE  SENSITIVITY  SUMMARY 


Parameter 

8 

a 

Fn 

GM 

Trim 

GZ(<f>) 

V 

A 

M 

N 

N 

L 

N 

B 

M 

N 

N 

M 

N 

C 

H 

N 

N 

H 

N 

r 

A 

N 

L 

N 

L 

N 

B 

N 

L 

N 

L 

N 

C 

L 

L 

N 

L 

N 

<t> 

A 

M 

H 

H 

X 

H 

B 

M 

H 

H 

H 

H 

C 

H 

H 

M 

H 

M 

coupled 

A 

M 

M 

X 

X 

X 

damping 

B 

L 

M 

H 

X 

H 

C 

X 

X 

X 

X 

H 

uncoupled 

A 

L 

H 

X 

X 

X 

damping 

B 

M 

H 

H 

X 

H 

C 

H 

H 

H 

X 

H 

This  reference  table  shows  clearly  that  the  sway  velocity  and  turning  rate  are 
relatively  insensitive  to  Froude  number,  metacentric  height  and  righting  moment, 
while  being  quite  sensitive  to  the  speed  loss  and  trim  parameters.  Roll  and  damping 
ratio  show  the  greatest  response  for  all  the  parameters  considered. 
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VI.   CONCLUSIONS  AND  RECOMMENDATIONS 

A.      CONCLUSIONS 

This  thesis  developed  a  mathematical  model  for  the  coupled  sway,  yaw  and  roll 
equations  of  motion  with  rudder  deflection  as  the  independent  variable.  Additionally, 
various  components  of  the  hydrodynamic  forces  and  moments,  e.g.,  rudder  to  hull 
interaction,  were  modelled  using  existing  research  formulation.  Deviation  of  this 
study  from  current  solution  methods  was  in  the  manner  of  evaluating  the  Taylor 
expansion  for  the  steady  state  Jacobian  matrix.  The  equilibrium  position  for  each 
rudder  angle  was  found  to  deviate  significantly  from  the  commonly  used  zero  state, 
even  for  rudder  angles  of  less  than  10°.  In  fact,  the  roll  angle  and  sway  velocity 
nominal  points  are  nearly  maximum  at  10°  of  rudder. 

Comparison  of  stability  results  for  the  coupled  lateral  motion  model  to  those 
of  models  representing  the  steering  equations  of  motion  decoupled  from  the  roll 
equations  of  motion  showed: 


1.  The  directional  stability  characteristics  signified  by  the  real  part  of  the 
complex  eigenvalues  are  not  affected  significantly  by  the  coupling  of  roll  into 
sway  and  yaw.  This  indicates  that  sway  and  yaw  can  be  treated  separately 
from  roll  when  studying  directional  stability. 

2.  Roll  stability  represented  by  the  magnitude  of  the  damping  ratio  is 
considerably  affected  by  the  coupling  of  sway  and  yaw  into  roll.  The  fact  that 
the  damping  ratio  is  underestimated  by  a  factor  of  two  by  the  uncoupled 
treatment  of  sway  and  yaw,  indicates  ship  design  by  this  means  is  a  significant 
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over-design    that   bears    more    attention,    especially   with    the    economic 
constraints  of  today's  Navy. 

3.  The  natural  frequency  is  not  significantly  affected  by  coupling.  This  indicates 
that  while  maneuvering  in  a  seaway  studies  require  the  coupling  of  the  lateral 
motions,  the  induced  roll  from  the  static  righting  moment  can  be  studied  from 
decoupled  systems. 

4.  Steering  stability  increases  with  increasing  rudder  deflection.  This  may 
indicate  that  if  a  ship  is  directionally  unstable  for  straight  line  motion,  it 
becomes  directionally  stable  during  a  turn  through  the  coupling  of  the  lateral 
motions. 

The  sensitivity  of  certain  design  parameters  to  the  coupled  system  was 

investigated  with  the  following  determinations: 


1.  The  parameters  tested  were  speed  loss  during  a  turn,  transverse  metacentric 
height,  Froude  number  of  initial  forward  speed,  trim  and  righting  moment.  Of 
these  parameters,  trim  proved  to  be  the  most  affected.  This  indicates  that 
during  all  phases  of  design,  operating  trim  must  be  considered  in  addition  to 
the  design  trim. 

2.  Sway  velocity  and  turning  rate  were  relatively  insensitive  for  all  parameters 
except  for  trim  and  speed  loss. 

3.  Roll  angle  and  damping  ratio  of  all  the  parameters  were  significantly  affected. 

4.  The  coupled  model  damping  ratio  showed  significantly  lower  sensitivity  than 
the  decoupled  model  for  small  rudder  angles  in  all  parameters. 


B.       RECOMMENDATIONS 

The  results  produced  in  this  thesis  were  as  a  consequence  of  experimentally 
reported  nonlinear  hydrodynamic  coefficients  and  design  parameters  of  a  SR-108 
container  ship  [Ref.  10].  Additional  studies  using  similar  data  are  required. 
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Utilizing  the  method  of  determining  equilibrium  points  for  each  rudder 
position,  design  a  control  law  for  an  autopilot  and  determine  the  effects  of  the 
coupled  lateral  motions  on  antopilot  dynamics. 

The  sharp  change  in  roll  angle  (<f>)  for  small  values  of  rudder  angle  (8)  suggests 
that  this  ship  would  be  an  excellent  candidate  for  rudder  roll  stabilizing  techniques. 

Finally,  as  indicated  in  Chapter  I,  Section  B,  use  the  coupled  sway,  yaw  and  roll 
equations  of  motion  and  build  a  mathematical  model  which  treats  wave  effects  in 
regular  and  irregular  seas.  Once  this  model  is  developed,  analyze  maneuvering  in 
waves  from  the  standpoint  of  biased  roll  oscillations.  Biased  roll  oscillations  are 
meant  as  roll  about  non-zero  heel  angle. 
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APPENDIX  A 


*********************************************************************** 


PROGRAM  COUPLED 

This  program  solves  the  coupled 
equations.  In  addition,  the  yaw- 
coupled  from  the  roll  equations  a 
results  with  the  fully  coupled  re 
into  several  subroutines,  both  or 
programs . 


yaw-sway-roll  ship  maneuvering 
sway  steering  equations  are  de- 
nd  solved  in  order  to  compare 
suits.   The  program  is  broken 
iginal  and  standard  commercial 


******************************** 
******************************** 


************************************* 
************************************* 

* 

*  The  INCLUDE  statement  is  used  to  pass  all  the  variables  between 

*  the  many  subroutines. 


* 

* 

** 

** 

* 
* 
* 
* 


INCLUDE  'HYDCOMN.for' 

COMMON /WORKSP/RWKSP 

REAL  RWKSP(5000) 

REAL*8  VDEL(300) ,V(300) ,R(300) ,RPHI(300) 

EXTERNAL  FCN,LSJAC,DNEQNJ 


Open  the  input  and  output  files 


0PEN(10,  FILE  ='COEF.INP' , statu 
0PEN(15,  FILE  ='VAR. OUT', status 
0PEN(20,  FILE  = 'NLIN. DAT' , statu 
0PEN(25,  FILE  ='EIGVAL1 .DAT' ,st 
0PEN(30,  FILE  =' EIGVAL2.DAT* ,st 
0PEN(26,  FILE  =' EIGVALL.DAT' ,st 
0PEN(27,  FILE  =' EIGVALp.DAT' ,st 
0PEN(35,  FILE  = 'ZETA. DAT ', statu 
0PEN(36,  FILE  =' OMEGA. DAT ', Stat 
0PEN(37,  FILE  ='ZR0LL.DAT' ,stat 
0PEN(38,  FILE  = 'RUDDER. DAT ' ,sta 
PI=4.0D0*DATAN(1.0D0) 
ERRREL=0.0001 
ITMAX=200 


s='old') 
='new' ) 
s='new' ) 
atus='new' ) 
atus='new' ) 
atus='new' ) 
atus='new' ) 
s='new* ) 
us='new' ) 
us='new' ) 
tus='new' ) 


This  subroutine  reads  the  hydrodynamic  derivative  and  coefficient 
values  from  COEF.INP  and  writes  pertinent  values  in  VAR.OUT. 

CALL  INPUT 
DEL=DEL*(PI/180.0D0) 

This  subroutine  sets  up  the  linear  EOM's  in  matrix  form: 

A  x  =  B 

CALL  LINEAR 

This  is  an  IMSL  MATH/LIBRARY  Subroutine  [Ref .  19] 
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Purpose:   Solve  a  real  general  system  of  linear  equations  with 


*  iterative  refinement 

* 

* 


Arguments : 


N     -Number  of  equations.   (Input) 

A     -N  by  N  matrix  containing  the  coefficients  of  the  linear 

system.   (Input) 
LDA   -Leading  dimension  of  A.   (Input) 

B     -Vector  of  length  N  containing  the  right-hand  side  of  the 
*  linear  system.   (Input) 

IPATH  -Path  indicator,  if  =  1  the  system  A*x=B  is  solved 
X     -Vector  of  length  N  containing  the  solution  to  the  linear 
system.   (Output) 


CALL  DLSARG(N,A,3,B,1,XGUESS) 

WRITE(15,*)'THE  SOLUTION  TO  THE  LINEAR  SYSTEM  OF  EQUATIONS  IS' 

WRITE(15,*)'  ' 

WRITE(15,*) 'V0=' ,XGUESS(1) 

WRITE(15,*) 'R0=' ,XGUESS(2) 

WRITE(15,*) 'PHIO=' ,XGUESS(3) 

WRITE(15,*)'  ' 

PHI0=XGUESS(3)*(180.0D0/PI) 

WRITE(15,*)'PHI0  IN  DEGREES  =*,PHIO 

WRITE(15,*)'  ' 

Set  the  rudder  limits  for  steady  turning 

DELMIN=0.01D0 
DELMAX=30.0D0 

This  Do-Loop  calculates  the  values  of  sway  velocity,  yaw  rate  and 
roll  angle  for  each  rudder  angle. 

DO  100  IDEL=1,IMAX 

WRITE  (*,*)  IDEL,IMAX 

RDEL=DELMIN+(DELMAX-DELMIN) *REAL( ( IDEL- 1 ) / ( IMAX- 1 ) ) 

DEL=RDEL*(PI/180.0D0) 

UND=1 .0D0-SLR*(DABS(DEL)) 

This  is  an  IMSL  MATH/LIBRARY  Subroutine  [Ref .  17] 

Purpose:   Solve  a  system  of  nonlinear  equations  using  the 

Levenberg-Marquardt  algorithm  with  a  user-supplied 

*  Jacobian. 

*  Arguments: 

* 

*  FCN    -  User-supplied  SUBROUTINE  to  evaluate  the  system  of 

*  equations  to  be  solved.   The  usage  is 

*  CALL  FCN  (X,  F,  N),  where 

*  X      -  The  point  at  which  the  functions  are 

*  evaluated.   (Input) 

*  X  should  not  be  changed  by  FCN. 
F      -  The  computed  function  values  at  the  point  X. 

(Output) 
N      -  Length  of  X,  F.   (Input) 

FCN  must  be  declared  EXTERNAL  in  the  calling  program. 
LSJAC   -  User-supplied  SUBROUTINE  to  evaluate  the  Jacobian  at  i 
point  X.   The  usage  is 
CALL  LSJAC  (N,  X,  FJAC) ,  where 
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N      -  Length  of  X.   (Input) 

X      -  The  point  at  which  the  function  is  evaluated. 
(Input) 

X  should  not  be  changed  by  LSJAC. 
FJAC   -  The  computed  N  by  N  Jacobian  at  the  point  X. 

(Output) 
LSJAC  must  be  declared  EXTERNAL  in  the  calling  program 
ERRREL  -  Stopping  criterion.   (Input) 

The  root  is  accepted  if  the  relative  error  between  two 
successive  approximations  to  this  root  is  less  than 
ERRREL. 
N      -  The  number  of  equations  to  be  solved  and  the  number 

of  unknowns.   (Input) 
ITMAX   -  The  maximum  allowable  number  of  iterations.   (Input) 
The  maximum  number  of  calls  to  FCN  is  ITMAX* (N+1). 
Suggested  value  =  200. 
XGUESS  -  A  vector  of  length  N.   (Input) 

XGUESS  contains  the  initial  estimate  of  the  root. 
X      -  A  vector  of  length  N.   (Output) 

X  contains  the  best  estimate  of  the  root  found  by 
DNEQNJ . 
FNORM   -  A  scalar  which  has  the  following  value, 

F(1)**2+.. .+F(N)**2  at  the  point  X.   (Output) 


CALL  DNEQNJ ( FCN , LSJAC , ERRREL , N , ITMAX , XGUESS , X , FNORM ) 
PHI=X(3)*180.0D0/PI 
WRITE ( 20,50 )RDEL,X(1),X( 2), PHI 
50   F0RMAT(4E15.5) 

Calculate  the  rudder  forces  and  moments:   Equations  2.6  and  2.7 

UP=UND* ( WP+TAU* ( (X ( 1 ) +XP*X ( 2 ) ) *  *2+CPV*X ( 1 ) +CPR*X ( 2 ) ) ) 

JP=(U*UP)/(NP*D) 

KT=KTA+(KTB*JP) 

UR=UP*EPSIL0N*((1.0D0+8.0D0*K*KT)/(PI*(JP**2)))**0.5 

VR=GAMMA*X(1)+CDR*X(2)+CDRRR*(X(2)**3)+CDRRV*X(1)*(X(2)**2) 

VD=(UR**2+VR**2)**0.5 

ALPHAR=DEL+DATAN(VR/UR) 

FLAM=(6.13D0*LAMBDA)/(LAMBDA+2.25D0) 

FNN=FLAM*(AR/(L**2) )*( (VD) **2) *DSIN(ALPHAR) 

YD=-(1.0D0+AH)*FNN*DC0S(DEL) 
ND=- (XR+AH*XH)*FNN*DCOS(DEL) 
KD= ( 1 . ODO+AH ) *ZR*FNN*DCOS ( DEL ) 
WRITE(38,*)ND,YD,KD 

Rudder  angle  in  degrees 

VDEL(IDEL)=DEL 

V(IDEL)=X(1) 

R(IDEL)=X(2) 

Rudder  angle  in  radians 

RPHI(IDEL)=X(3) 
XGUESS(1)=X(1) 
XGUESS(2)=X(2) 
XGUESS(3)=X(3) 
100   CONTINUE 
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* 


* 


This  subroutine  solves  the  eigenvalue  problem  |A  -  lambda*B|=  0 
for  the  fully  coupled  yaw-sway-roll  equations,  using  the  values 
for  v,  r,  and  phi  obtained  above. 

CALL  EIGTH1(VDEL,V,R,RPHI) 

This  subroutine  solves  the  eigenvalue  problem  |As  -  lambda*Bs|=  0 
for  the  steering  yaw-sway  equations. 

CALL  EIGTH2(VDEL,V,R) 

This  subroutine  solves  the  eigenvalue  problem  |Ar  -  lambda*Br|=  0 
for  the  roll  equations. 

CALL  EIGTH3(VDEL,Rphi) 

STOP 

END 


SUBROUTINE  INPUT 
INCLUDE  'HYDCOMN.for' 

*  The  numbers  in  the  variables  such  as  N122  stand  for 

*  1  implies  v 

*  2  implies  r 

*  3  implies  phi 

*  4  implies  delta 


READ(10,*)N0,N1,N2,N3,N111,N122,N222,N112,N113,N133,N223, 
&  N233 , KO , K1 , K2 , K3 , K1 1 1 , K1 22 , K222 ,  K1 1 2 ,  K1 1 3 ,  K1 33 , K223 , K233 , 
&  Y0,Y1,Y2,Y3,Y111,Y122,Y222,Y112JY113,Y133,Y223,Y233,L,W,XG, 
&  AH, AR, LAMBDA, XR ,XH,ZR ,BM,WP,TAU,XP,CPV ,CPR ,EPSILON,K,KTA,KTB, 
&  D,CDR,CDRRR,CDRRV, GAMMA 


* 

*   Calculate  the  nondimensional  mass 


M=W/(0.5D0*(L**3)) 

WRITE (*,*) 'ENTER  THE  GM  AND  THE  SPEED  LOSS  RATIO' 

READ(*,*)GM,SLR 

WRITE(15,*)'GM=' ,GM 

WRITE(15,*)'  ' 

WRITE(15,*)'THE  SPEED  LOSS  RATIO  IS=' ,SLR 

Adjusts  the  z  coordinate  of  CG  based  on  the  GM 

ZGM=DABS(0.3D0-GM)/L 

ZG=0. 00078 -ZGM 

WRITE(15,*)'THE  SHIP  DISPLACEMENT  (IN  CUBIC  METERS)  IS  W=*,W 

WRITE(15,*)'  ' 

WRITE(15,*)'THE  SHIP  NONDIMENSIONAL  MASS  IS  M=',M 

WRITE(15,*)'  ' 

XG=XG/L 

WRITE(15,*)'L=',L 

WRITE(15,*)'XG=',XG 

WRITE(15,*)'ZG=',ZG 

WRITE(15,*)'  ' 

WRITE (*,*) 'ENTER  THE  INITIAL  RUDDER  ANGLE  IN  DEGREES' 

READ(*,*)DEL 

WRITE(15,*)'THE  INITIAL  RUDDER  ANGLE  DEL=',DEL 

WRITE (*,*) 'ENTER  THE  FROUDE  NUMBER  AND  THE  IMAX' 

READ(*,*)FN,IMAX 
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WRITE(15,*)'THE  FROUDE  NUMBER  IS  FN=',FN 

Calculate  the  initial  forward  speed 

U=FN*((9.8D0*L)**0.5D0) 

WRITE(15,*)'THE  FORWARD  VELOCITY  IS  U=',U 

WRITE(15,*)'  ' 

WRITE (*,*)' ENTER  THE  PROPELLER  SPEED  IN  RPMS  FOR  THE  FN' 

READ(*,*)NP 

WRITE(15,*)'THE  PROPELLER  SPEED  IN  RPMS  IS  NP=',NP 

Change  rpm's  to  radians  per  sec 

NP=NP*(PI/30.0D0) 


The  weight  component  is  nondimensionalized  such  that  WG*GM  is 
*   dimensionless 


* 


WG=(W*9.81D0)/(0.5D0*(L**3)*(U**2)) 

WRITE(15,*)'THE  NONDIM  COMPONENT  OF  DISPLACEMENT  IS  WG=',WG 

RETURN 

END 

SUBROUTINE  LINEAR 


INCLUDE  'HYDCOMN.for' 

INTEGER  IROW,  ICOL 

DO  100  IR0W=1,N 

B(IR0W)=0.0D0 

DO  90  IC0L=1,N 

A( IROW, ICOL )=0.0D0 

90 

CONTINUE 

* 

100 

CONTINUE 

* 
* 

I 

Equation  3.2 

A(1,1)=N1 

A(1,2)=N2-(M*XG*UND) 

A(1,3)=N3 

A(2,1)=Y1 

A(2,2)=Y2-(M*UND) 

A(2,3)=Y3 

A(3,1)=K1 

A(3,2)=K2+(M*ZG*UND) 

A(3,3)=K3-(WG*GM) 

FNL=((6.13D0*LAMBDA)/(LAMBDA+2.25D0))*(AR/L**2) 

N4=-(XR+(AH*XH))*FNL 

Y4=-(1.0D0+AH)*FNL 

K4=(1.0D0+AH)*ZR*FNL 

B(1)=-(N4*DEL)-N0 

B(2)=-(Y4*DEL)-Y0 

B(3)=-(K4*DEL)-K0 

RETURN 

END 

SUBROUTINE  FCN(Z,Q,NDIM) 
INCLUDE  'HYDCOMN.for* 
INTEGER  NDIM 
DIMENSION  Z(NDIM),Q(NDIM) 
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Equation  2.7 

UP=UND*(WP+TAU*( (Z(1 )+XP*Z(2) )**2+CPV*Z(1 )+CPR*Z(2) ) ) 

JP=(U*UP)/(NP*D) 

KT=KTA+(KTB*JP) 

UR=UP*EPSIL0N*((1.0D0+8.0D0*K*KT)/(PI*(JP**2)))**0.5 

VR=GAMMA*Z(1)+CDR*Z(2)+CDRRR*(Z(2)**3)+CDRRV*Z(1)*(Z(2)**2) 

VD=(UR**2+VR**2)**0.5 

ALPHAR=DEL+DATAN(VR/UR) 

FLAM=( 6. 13D0* LAMBDA) / (LAMBDA+2.25D0) 

FNN=FLAM*(AR/(L**2))*((VD)**2)*DSIN(ALPHAR) 

Equation  2.6 

Y4=-(1.0D0+AH)*FNN*DC0S(DEL) 
N4=- (XR+AH*XH)*FNN*DCOS(DEL) 
K4=(1 .ODO+AH)*ZR*FNN*DCOS(DEL) 

Equation  2.10b 

C1=N1*Z(1)+N111*(Z(1)**3)+N122*Z(1)*(Z(2)**2) 
C2=N222*(Z(2)**3)+N112*(Z(1)**2)*Z(2)+N3*Z(3) 

C3=N113*(Z(1)**2)*Z(3)+N133*Z(1)*(Z(3)**2)+(N2-(M*XG*UND))*Z(2) 
C4=N233*Z(2)*(Z(3)**2)+N4+N0+N223*(Z(2)**2)*Z(3) 

Q(1)=C1+C2+C3+C4 

Equation  2.10c 

D1=Y1*Z(1)+Y111*(Z(1)**3)+Y122*Z(1)*(Z(2)**2)+(Y2-(M*UND))*Z(2) 
D2=Y222*(Z(2)**3)+Y112*(Z(1)**2)*Z(2)+Y3*Z(3) 

D3=Y113*(Z(1)**2)*Z(3)+Y133*Z(1)*(Z(3)**2)+Y223*(Z(2)**2)*Z(3) 
D4=Y233*Z(2)*(Z(3)**2)+Y4+Y0 

Q(2)=D1+D2+D3+D4 

Equation  2.9 

GZ=GM*DSIN(Z(3))+0.5*BM*(DTAN(Z(3)))**2*DSIN(Z(3)) 
GZ=GM*DSIN(Z(3)) 

Equation  2.10d 

E1=K1*Z(1)+K111*(Z(1)**3)+K122*Z(1)*(Z(2)**2)+(K2+M*ZG*UND)*Z(2) 
E2=K222*(Z(2)**3)+K112*(Z(1)**2)*Z(2)+K3*Z(3) 

E3=K113*(Z(1)**2)*Z(3)+K133*Z(1)*(Z(3)**2)+K223*(Z(2)**2)*Z(3) 
E4=K233*Z(2)*(Z(3)**2)+K4+K0-(WG*GZ) 

Q(3)=E1+E2+E3+E4 

RETURN 
END 

SUBROUTINE  LSJAC(NDIM,Z,QJAC) 

INCLUDE  'HYDCOMN.for' 

INTEGER  NDIM 

DIMENSION  Z(NDIM) ,QJAC(NDIM,NDIM) 

Equation  2.7 
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UP=UND*(WP+TAU*( (Z(1 )+XP*Z(2) )**2+CPV*Z(1 )+CPR*Z(2) ) ) 

V9=Z(1)*0.99D0 

R9=Z(2)*0.99D0 

UPV9=UND*(WP+TAU*((V9+XP*Z(2))**2+CPV*V9+CPR*Z(2))) 

UPR9=UND*(WP+TAU*((Z(1)+XP*R9)**2+CPV*Z(1)+CPR*R9)) 

JP=(U*UP)/(NP*D) 

JPV9=(U*UPV9)/(NP*D) 

JPR9=(U*UPR9)/(NP*D) 

KT=KTA+(KTB*JP) 

KTV9=KTA+(KTB*JPV9) 

KTR9=KTA+ ( KTB* J PR 9 ) 

UR=UP*EPSIL0N*((1.0D0+8.0D0*K*KT)/(PI*(JP**2)))**0.5 

URV9=UPV9*EPSIL0N*((1.0D0+8.0D0*K*KTV9)/(PI*(JPV9**2)))**0.5 

URR9=UPR9*EPSIL0N*( (1 .0D0+8.0D0*K*KTR9) / (PI* (JPR9**2) ) )**0.5 

VR=GAMMA*Z(1)+CDR*Z(2)+CDRRR*(Z(2)**3)+CDRRV*Z(1)*(Z(2)**2) 

VRV9=GAMMA*V9+CDR*Z(2)+CDRRR*(Z(2)**3)+CDRRV*V9*(Z(2)**2) 

VRR9=GAMMA*Z(1)+CDR*R9+CDRRR*(R9**3)+CDRRV*Z(1)*(R9**2) 

VD=(UR**2+VR**2)**0.5 

VDV9=(URV9**2+VRV9**2)**0.5 

VDR9=(URR9**2+VRR9**2)**0.5 

ALPHAR=DEL+DATAN(VR/UR) 

ALPHARV=DEL+DATAN(VRV9/URV9) 

ALPHARR=DEL+DATAN(VRR9/URR9) 

FLAM= (6 . 1 3D0* LAMBDA) / ( LAMBDA+2 . 25D0 ) 

FNN=FLAM*(AR/(L**2) )*( (VD) **2) *DSIN(ALPHAR) 

FNNV9=FLAM*(AR/(L**2))*((VDV9)**2)*DSIN(ALPHARV) 

FNNR9=FLAM*(AR/(L**2))*((VDR9)**2)*DSIN(ALPHARR) 

Y4=-(1.0D0+AH)*FNN*DC0S(DEL) 

Y4V9=-(1.0D0+AH)*FNNV9*DC0S(DEL) 

Y4R9=-(1.0D0+AH)*FNNR9*DC0S(DEL) 

*  Equation  3.4 

* 

PY4V=(Y4-Y4V9)/(Z(1)-V9) 

PY4R=(Y4-Y4R9)/(Z(2)-R9) 

N4=-(XR+AH*XH)*FNN*DC0S(DEL) 

N4V9=-(XR+AH*XH)*FNNV9*DC0S(DEL) 

N4R9=-(XR+AH*XH)*FNNR9*DC0S(DEL) 

PN4V=(N4-N4V9)/(Z(1)-V9) 

PY4R=(Y4-Y4R9)/(Z(2)-R9) 

N4=-(XR+AH*XH)*FNN*DC0S(DEL) 

*  Equation  3.4 

* 

N4V9=-(XR+AH*XH)*FNNV9*DC0S(DEL) 

N4R9=-(XR+AH*XH)*FNNR9*DC0S(DEL) 

PN4V=(N4-N4V9) / (Z( 1 ) -V9) 

PN4R=(N4-N4R9)/(Z(2)-R9) 

K4=(1 .ODO+AH)*ZR*FNN*DCOS(DEL) 

K4V9=(1.0D0+AH)*ZR*FNNV9*DC0S(DEL) 

K4R9= ( 1 . ODO+AH ) *ZR*FNNR9*DC0S ( DEL ) 

* 

*  Equation  3.4 


PK4V=(K4-K4V9) / (Z(1 ) -V9) 

PK4R=(K4-K4R9)/(Z(2)-R9) 

* 

*  Equation  3.3b 


C11=N1+3*N111*(Z(1)**2)+N122*(Z(2)**2)+2.0D0*N112*Z(1)*Z(2) 
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C12=2.0D0*N113*Z(1)*Z(3)+N133*(Z(3)**2) 

QJAC(1,1)=C11+C12+PN4V 

Equation  3.3c 

C21=2.0D0*N122*Z(1)*Z(2)+(N2-M*XG*UND)+3.0D0*N222*(Z(2)**2) 
C22=N112*(Z(1)**2)+2.0D0*N223*Z(2)*Z(3)+N233*(Z(3)**2) 

QJAC(1 ,2)=C21+C22+PN4R 

Equation  3.3d 

C31=N3+N113*(Z(1)**2)+2.0D0*N133*Z(1)*Z(3)+N223*(Z(2)**2) 
C32=2.0D0*N233*Z(2)*Z(3) 

QJAC(1,3)=C31+C32 

Equation  3.3e 

D11=Y1+3.0D0*Y111*(Z(1)**2)+Y122*(Z(2)**2)+2.0D0*Y112*Z(1)*Z(2) 
D12=2.0D0*Y113*Z(1)*Z(3)+Y133*(Z(3)**2) 

QJAC(2,1)=D11+D12+PY4V 

Equation  3.3f 

D21=2.0D0*Y122*Z(1)*Z(2)+(Y2-M*UND)+3.0D0*Y222*(Z(2)**2) 
D22=2.0D0*Y223*Z(2)*Z(3)+Y233*(Z(3)**2)+Y112*(Z(1)**2) 

QJAC(2,2)=D21+D22+PY4R 

Equation  3.3g 

D31=Y3+Y113*(Z(1)**2)+2.0D0*Y133*Z(1)*Z(3)+Y223*(Z(2)**2) 
D32=2.0D0*Y233*Z(2)*Z(3) 

QJAC(2,3)=D31+D32 

Equation  3.3h 

E11=K1+3.0D0*K111*(Z(1)**2)+K122*(Z(2)**2)+2.0D0*K112*Z(1)*Z(2) 
E12=2.0D0*K113*Z(1)*Z(3)+K133*(Z(3)**2) 

QJAC(3,1)=E11+E12+PK4V 

Equation  3.3i 

E21=2.0D0*K122*Z(1)*Z(2)+(K2+M*ZG*UND)+3.0D0*K222*(Z(2)**2) 
E22=2.0D0*K223*Z(2)*Z(3)+K233*(Z(3)**2)+K112*(Z(1)**2) 

QJAC(3,2)=E21+E22+PK4R 

Equation  3.3k 


c      DGZ1=GM*DC0S(Z(3))+0.5D0*BM*(DTAN(Z(3)))**2*DC0S(Z(3)) 
c      DGZ2=BM*(DTAN(Z(3))*((1.0D0/DC0S(Z(3)))**2)*DSIN(Z(3))) 
c      DGZ=DGZ1+DGZ2 

DGZ=GM*DC0S(Z(3)) 


Equation  3.3] 
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E31=K3+K113*(Z(1)**2)+2.0D0*K133*Z(1)*Z(3)+K223*(Z(2)**2) 
E32=2.0D0*K233*Z(2)*Z(3)-WG*DGZ 

QJAC(3,3)=E31+E32 

RETURN 
END 

SUBROUTINE  EIGTH1 (VDEL,V,R,RPHI) 

INCLUDE  'HYDCOMN.for' 
INTEGER  I,J,IERR 

REAL*8  AE(4,4) ,BE(4,4) ,BETA(4) ,VDEL(300) ,V(300) ,R(300) ,RPHI(300; 
REAL*8  P(300) ,EVALR(4) ,EVALI(4) ,EVAL(4) ,ALFR(4) ,ALFI(4) ,Z(4,4) 
EXTERNAL  DMACH,RGG 
DO  500  I=1,IMAX 
P(I)=O.ODO 
SET  UP  THE  B  MATRIX  AND  RESET  EACH  I 

UP=UND* (WP+TAU* ( ( V ( I ) +XP*R ( I ) ) **2+CPV*V ( I ) +CPR*R ( I ) ) ) 

V9=V(I)*0.99D0 

R9=R ( I ) *0 . 99D0 

UPV9=UND* (WP+TAU* ( ( V9+XP*R ( I ) )  **2+CPV*V9+CPR*R ( I ) ) ) 

UPR9=UND* (WP+TAU* ( (V( I )+XP*R9) **2+CPV*V( I ) +CPR*R9) ) 

JP=(U*UP)/(NP*D) 

JPV9=(U*UPV9)/(NP*D) 

JPR9=(U*UPR9)/(NP*D) 

KT=KTA+(KTB*JP) 

KTV9=KTA+(KTB*JPV9) 

KTR9=KTA+ ( KTB* JPR9 ) 

UR=UP*EPSIL0N*((1.0D0+8.0D0*K*KT)/(PI*(JP**2)))**0.5 

URV9=UPV9*EPSIL0N*((1.0D0+8.0D0*K*KTV9)/(PI*(JPV9**2)))**0.5 

URR9=UPR9*EPSIL0N*((1.0D0+8.0D0*K*KTR9)/(PI*(JPR9**2)))**0.5 

VR=GAMMA* V ( I ) +CDR  *R ( I ) +CDRRR  * ( R ( I ) *  *3 ) +CDRRV*V (I)*(R(I)**2) 

VRV9=GAMMA*V9+CDR*R ( I ) +CDRRR* ( R ( I ) *  *3 ) +CDRRV*V9* ( R ( I ) **2 ) 

VRR9=GAMMA*V(I)+CDR*R9+CDRRR*(R9**3)+CDRRV*V(I)*(R9**2) 

VD=(UR**2+VR**2)**0.5 

VDV9=(URV9**2+VRV9**2)**0.5 

VDR9=(URR9**2+VRR9**2)**0.5 

ALPHAR=VDEL(I)+DATAN(VR/UR) 

ALPHARV=VDEL(I)+DATAN(VRV9/URV9) 

ALPHARR=VDEL ( I ) +DATAN ( VRR9/URR9 ) 

FLAM=(6.13D0*LAMBDA)/(LAMBDA+2.25D0) 

FNN=FLAM*(AR/(L**2) )*( (VD) **2) *DSIN (ALPHAR) 

FNNV9=FLAM*(AR/(L**2))*((VDV9)**2)*DSIN(ALPHARV) 

FNNR9=FLAM*(AR/(L**2))*((VDR9)**2)*DSIN(ALPHARR) 

Y4=-(1.0D0+AH)*FNN*DC0S(VDEL(I)) 

Y4V9=-(1.0D0+AH)*FNNV9*DC0S(VDEL(I)) 

Y4R9=-(1.0D0+AH)*FNNR9*DC0S(VDEL(I)) 

PY4V=(Y4-Y4V9) / (V( I ) -V9) 

PY4R=(Y4-Y4R9)/(R(I)-R9) 

N4=-(XR+AH*XH)*FNN*DC0S(VDEL(I)) 

N4V9=-(XR+AH*XH)*FNNV9*DC0S(VDEL(I)) 

N4R9=-(XR+AH*XH)*FNNR9*DC0S(VDEL(I)) 

PN4V=(N4-N4V9)/(V(I)-V9) 

PN4R=(N4-N4R9)/(R(I)-R9) 

K4= ( 1 . ODO+AH ) *ZR*FNN*DCOS ( VDEL ( I ) ) 

K4V9=(1.0D0+AH)*ZR*FNNV9*DC0S(VDEL(I)) 

K4R9=(1.0D0+AH)*ZR*FNNR9*DC0S(VDEL(I)) 

PK4V=(K4-K4V9) / (V( I ) -V9) 
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PK4R=(K4-K4R9)/(R(I)-R9) 

Calculate  the  fundamental  circular  frequency 

WN=((WG*GM)/(0.000021D0))**0.5 

IF(FN.LE.0.1)THEN 

KAPPA=0.1D0 

ELSE  IF(FN.GE.0.2)THEN 

KAPPA=0.2D0 

ELSE 

KAPPA=FN 

ENDIF 

KPHID0T=-0. 000021 DO*KAPPA*WN 


Equation  4.1b 


BE(1 
BE(1 
BE(1 
BE(1 
BE(2 
BE(2 
BE(2 
BE(2 
BE(3 
BE(3 
BE(3 
BE(3 
BE(4 
BE(4 
BE(4 
BE(4 


=0.00035245D0 

=0.000875D0 

=0.0D0 

=-0.000213D0 

=0.014969D0 

=0.00035245D0 

=-0.000221D0 

=0.0D0 

=-0.000221D0 

=0.0D0 

=0.000021  DO 

=-KPHID0T 

=0.0D0 

=0.0D0 

=0.0D0 

=1.0D0 


Equation  4.1c 


SET  UP  THE  A  MATRIX 


C11=N1+3*N111*(V(I)**2)+N122*(R(I)**2)+2.0D0*N112*V(I)*R(I) 
C12=2.0D0*N113*V(I)*RPHI(I)+N133*(RPHI(I)**2) 

AE(1,1)=C11+C12+PN4V 

C21=2.0D0*N122*V(I)*R(I)+(N2-M*XG*UND)+3.0D0*N222*(R(I)**2) 
C22=N1 1 2* ( V ( I ) **2) +2 . 0D0*N223*R ( I ) *RPHI ( I ) +N233* (RPHI ( I ) **2 ) 

AE(1,2)=C21+C22+PN4R 
AE(1,3)=0.0D0 

C41=N3+N113*(V(I)**2)+2.0D0*N133*V(I)*RPHI(I)+N223*(R(I)**2) 
C42=2 .0D0*N233*R ( I ) *RPHI ( I ) 

AE(1,4)=C41+C42 

D11=Y1+3.0D0*Y111*(V(I)**2)+Y122*(R(I)**2)+2.0D0*Y112*V(I)*R(I) 
D12=2.0D0*Y113*V(I)*RPHI(I)+Y133*(RPHI(I)**2) 

AE(2,1)=D11+D12+PY4V 

D21 =2 . 0D0*Y1 22*V ( I ) *R ( I ) + ( Y2-M*UND) +3 . 0D0*Y222* (R ( I ) **2) 
D22=2 .0D0*Y223*R ( I ) *RPHI ( I )+Y233* (RPHI ( I ) **2) +Y1 12* (V( I ) **2) 

AE(2,2)=D21+D22+PY4R 
AE(2,3)=0.0D0 
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D41=Y3+Y113*(V(I)**2)+2.0D0*Y133*V(I)*RPHI(I)+Y223*(R(I)**2) 
D42=2 . 0D0*Y233*R ( I ) *RPHI ( I ) 

AE(2,4)=D41+D42 

E11=K1+3.0D0*K111*(V(I)**2)+K122*(R(I)**2)+2.0D0*K112*V(I)*R(I) 
E12=2.0D0*K113*V(I)*RPHI(I)+K133*(RPHI(I)**2) 

AE(3,1)=E11+E12+PK4V 

E21=2.0D0*K122*V(I)*R(I)+(K2+M*ZG*UND)+3.0D0*K222*(R(I)**2) 
E22=2.0D0*K223*R(I)*RPHI(I)+K233*(RPHI(I)**2)+K112*(V(I)**2) 

AE(3,2)=E21+E22+PK4R 
AE(3,3)=0.0D0 

DGZ1=GM*DC0S(RPHI(I))+0.5D0*BM*(DTAN(RPHI(I)))**2*DC0S(RPHI(I)) 

DGZ2=BM* ( DTAN ( RPHI ( I ) ) * ( ( 1 . ODO / DCOS ( RPH I ( I ) ) ) *  *2 ) *DS IN ( RPH I ( I ) ) ) 

DGZ=DGZ1+DGZ2 

DGZ=GM*DCOS(RPHI(I)) 
E41=K3+K113*(V(I)**2)+2.0D0*K133*V(I)*RPHI(I)+K223*(R(I)**2) 
E42=2 .0D0*K233*R ( I ) *RPHI ( I ) -WG*DGZ 

AE(3,4)=E41+E42 

AE(4,1)=0.0D0 

AE(4,2)=0.0D0 

AE(4,3)=1.0D0 

AE(4,4)=0.0D0 

This  is  a  subroutine  that  utilizes  the  EISPACK  LIBRARY  [Ref.  20] 

Purpose:   Calls  the  recommended  sequence  of  subroutines  from  the 
Eigensystem  subroutine  package  (EISPACK)  to  find  the 
eigenvalues  of  the  real  generalized  eigenproblem: 
AX  =  ( LAMBDA) *BX 


Arguments: 


NM 


N 

A 

B 

MATZ 

ALFR 

ALFI 

BETA 

Z 
IERR 


Row  dimension 
declared  in  t 
(Input) 
The  order  of 
Contains  a  re 
Contains  a  re 
Path  indicato 
Real  parts  of 
•Imaginary  par 
(Output) 
•Contains  the 
given  by  the 
•For  eigenvect 
■Error  complet 


of  the  two-dimensional  array  parameters  as 
he  calling  program  Dimension  statement 

the  matrices  A  and  B.    (Input) 

al  general  matrix   (Input) 

al  general  matrix   (Input) 

r.   If=  0,  eigenvalues  are  found.   (Input) 

the  numerators  of  the  eigenvalues.  (Output) 
ts  of  the  numerators  of  the  eigenvalues. 

denominators  of  the  eigenvalues,  which  are 
ratios:   (ALFR+I*ALFI) /BETA.   (Output) 
or  usage  if  MATZ  not  equal  to  zero.  (Output) 
ion  code.   (Output) 


CALL  RGG(4, 4, AE, BE, ALFR, ALFI, BETA, 0,Z, IERR) 
DO  600  J=1,4 
IF(BETA(J) .NE.0.0D0)THEN 

EVALR (J ) =ALFR ( J ) /BETA( J ) 

EVALI (J ) =ALFI (J ) /BETA( J ) 
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ELSE 

EVALR(J)=DMACH(2) 
EVALI(J)=DMACH(2) 
ENDIF 
600   CONTINUE 

WRITE(25,25)EVALR(1),EVALI(1) , EVALR (2) , EVALI (2) 
WRITE(30,25)EVALR(3) ,EVALI(3) ,EVALR(4) ,EVALI(4) 
* 

*  Calculate  the  damping  ratio  and  the  natural  frequency 

* 

Z1=- EVALR (1)/(DSQRT( EVALI ( 1 ) **2+EVALR ( 1 ) **2) ) 
Z2=-EVALR(2)/(DSQRT(EVALI(2)**2+EVALR(2)**2)) 
Z3=-EVALR(3)/(DSQRT(EVALI(3)**2+EVALR(3)**2)) 
Z4=-EVALR(4)/(DSQRT(EVALI(4)**2+EVALR(4)**2)) 

WN1 =DSQRT( EVALI ( 1 ) **2+EVALR ( 1 ) **2) 

WN2=DSQRT(EVALI(2)**2+EVALR(2)**2) 

WN3=DSQRT(EVALI(3)**2+EVALR(3)**2) 

WN4=DSQRT(EVALI(4)**2+EVALR(4)**2) 

WRITE(35,25)Z1,Z2,Z3,Z4 

WRITE(36,25)WN1 ,WN2,WN3,WN4 
25   F0RMAT(4E15.5) 
500  CONTINUE 

RETURN 

END 

* 

SUBROUTINE  EIGTH2(VDEL,V,R) 

INCLUDE  'HYDCOMN.for' 

INTEGER  I,J,IERR 

REAL*8  AE(2,2) ,BE(2,2) ,BETA(2) ,VDEL(300) ,V(300) ,R(300) 

REAL*8  EVALR(2) ,EVALI(2) ,EVAL(2) ,ALFR(2) ,ALFI(2) ,Z(2,2) 

EXTERNAL  DMACH.RGG 

DO  500  1=1 ,imax 

*  SET  UP  THE  B  MATRIX  AND  RESET  EACH  I 

UP=UND* (WP+TAU* ( ( V ( I ) +XP*R ( I ) ) **2+CPV*V ( I ) +CPR*R ( I ) ) ) 

V9=V(I)*0.99D0 

R9=R ( I ) *0 .99D0 

UPV9=UND* (WP+TAU* ( ( V9+XP*R ( I ) ) **2+CPV*V9+CPR*R ( I ) ) ) 

UPR9=UND*(WP+TAU*((V(I)+XP*R9)**2+CPV*V(I)+CPR*R9)) 

JP=(U*UP)/(NP*D) 

JPV9=(U*UPV9)/(NP*D) 

JPR9=(U*UPR9)/(NP*D) 

KT=KTA+(KTB*JP) 

KTV9=KTA+(KTB*JPV9) 

KTR9=KTA+ (KTB*JPR9) 

UR=UP*EPSIL0N*((1.0D0+8.0D0*K*KT)/(PI*(JP**2)))**0.5 

URV9=UPV9*EPSIL0N*((1.0D0+8.0D0*K*KTV9)/(PI*(JPV9**2)))**0.5 

URR9=UPR9*EPSIL0N*( (1 .0D0+8.0D0*K*KTR9) / (PI* (JPR9**2) ) )**0.5 

VR=GAMMA*V ( I ) +CDR  *R ( I ) +CDRRR  * ( R ( I ) *  *3 ) +CDRRV*V (I)*(R(I)**2) 

VRV9=GAMMA*V9+CDR*R ( I ) +CDRRR* (R ( I ) **3) +CDRRV*V9* (R ( I ) **2 ) 

VRR9=GAMMA*V( I )+CDR*R9+CDRRR* (R9**3) +CDRRV*V( I ) * (R9**2) 

VD=(UR**2+VR**2)**0.5 

VDV9= (URV9**2+VRV9**2) **0 . 5 

VDR9=(URR9**2+VRR9**2)**0.5 

ALPHAR=VDEL(I)+DATAN(VR/UR) 

ALPHARV=VDEL ( I ) +DATAN ( VRV9 / URV9 ) 

ALPHARR=VDEL(I)+DATAN(VRR9/URR9) 

FLAM=(6.13D0*LAMBDA)/(LAMBDA+2.25D0) 

FNN=FLAM*(AR/(L**2) )*( (VD)**2)*DSIN(ALPHAR) 

FNNV9=FLAM*(AR/(L**2))*((VDV9)**2)*DSIN(ALPHARV) 

FNNR9=FLAM*(AR/(L**2))*((VDR9)**2)*DSIN(ALPHARR) 
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Y4=-(1 .ODO+AH)*FNN*DCOS(VDEL(I)) 

Y4V9=-(1 .0D0+AH)*FNNV9*DC0S(VDEL(I)) 

Y4R9=-(1 .0D0+AH)*FNNR9*DC0S(VDEL(I)) 

PY4V=(Y4-Y4V9)/(V(I)-V9) 

PY4R=(Y4-Y4R9)/(R(I)-R9) 

N4=-(XR+AH*XH)*FNN*DC0S(VDEL(I)) 

N4V9=-(XR+AH*XH)*FNNV9*DC0S(VDEL(I)) 

N4R9=-(XR+AH*XH)*FNNR9*DC0S(VDEL(I)) 

PN4V=(N4-N4V9)/(V(I)-V9) 

PN4R=(N4-N4R9)/(R(I)-R9) 

WN= ( (WG*GM)/ (0.000021  DO ))**0. 5 

Equation  4.4b 


BE(1 ,1)=0.00035245D0 
BE(1,2)=0.000875D0 
BE(2,1)=0.014969D0 
BE(2,2)=0.00035245D0 

r 

r   Equation  4.4c-g 

r 

C11=N1+3*N111*(V(I)**2)+N122*(R(I)**2)+2.0D0*N112*V(I)*R(I) 

r 

AE(1,1)=C11+PN4V 

r 

C21=2.0D0*N122*V(I)*R(I)+(N2-M*XG*UND)+3.0D0*N222*(R(I)**2) 
C22=N112*(V(I)**2)+PN4R 

r 

AE(1,2)=C21+C22 
D11=Y1+3.0D0*Y111*(V(I)**2)+Y122*(R(I)**2)+2.0D0*Y112*V(I)*R(I) 

r 

AE(2,1)=D11+PY4V 

r 

D21=2.0D0*Y122*V(I)*R(I)+(Y2-M*UND)+3.0D0*Y222*(R(I)**2) 
D22=Y112*(V(I)**2)+PY4R 

r 

AE(2,2)=D21+D22 

r 

CALL  RGG(2,2,AE)BEJALFR,ALFI,BETA,0,Z,IERR) 
DO  600  J=1,2 
IF(BETA(J) .NE.0.0D0)THEN 
EVALR(J)=ALFR(J)/BETA(J) 
EVALI(J)=ALFI(J)/BETA(J) 
ELSE 

EVALR(J)=DMACH(2) 
EVALI(J)=DMACH(2) 
ENDIF 
600  CONTINUE 

WRITE(26,25)EVALR(1),EVALI(1),EVALR(2),EVALI(2) 
25   F0RMAT(4E15.5) 
500  CONTINUE 
RETURN 
END 

SUBROUTINE  EIGTH3(VDEL ,RPHI ) 

INCLUDE  'HYDCOMN.for' 
INTEGER  I,J,IERR 
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REAL*8  AE(2,2) ,BE(2,2) ,BETA(2) ,VDEL(300) ,RPHI(300) 
REAL*8  P(300) ,EVALR(2) ,EVALI(2) ,EVAL(2) ,ALFR(2) ,ALFI(2) ,Z(2,2) 
EXTERNAL  DMACH,RGG 
DO  500  1=1 ,imax 
P(I)=O.ODO 
SET  UP  THE  B  MATRIX  AND  RESET  EACH  I 
WN=( (WG*GM)/(0.000021D0) )**0.5 
IF(FN.LE.0.1)THEN 
KAPPA=0.1D0 

ELSE  IF(FN.GE.0.2)THEN 
KAPPA=0 . 2D0 
ELSE 

KAPPA=FN 
ENDIF 
KPHID0T= - 0 . 000021 DO*KAPPA*WN 

Equation  4.7b 

BE(1,1)=0.000021D0 
BE(1,2)=-KPHID0T 
BE(2,1)=0.0d0 
BE(2,2)=1.0D0 

Equation  4.7c-e 


AE(1,1)=O.OdO 
c      DGZ1 =GM*DCOS ( RPHI ( I ) ) +0 . 5D0*BM* ( DTAN ( RPHI ( I ) ) ) *  *2*DC0S ( RPHI ( I ) ) 
c      DGZ2=BM* (DTAN ( RPHI ( I ) ) * ( ( 1 . ODO/DCOS ( RPHI ( I ) ) ) **2 ) *DSIN (RPHI ( I ) ) ) 
c      DGZ=DGZ1+DGZ2 

DGZ=GM*DCOS(RPHI(I)) 
AE(1,2)=K3-WG*DGZ 
AE(2,1)=1.0D0 
AE(2,2)=0.0D0 

CALL  RGG(2,2,AE,BE,ALFR,ALFI,BETA,0,Z,IERR) 
DO  600  J=1 ,2 
IF(BETA(J) ,NE.O.ODO)THEN 
EVALR (J ) =ALFR ( J ) /BETA( J ) 
EVALI(J)=ALFI(J)/BETA(J) 
ELSE 

EVALR(J)=DMACH(2) 
EVALI(J)=DMACH(2) 
ENDIF 
600  CONTINUE 

WRITE(27,25)EVALR(1),EVALI(1),EVALR(2),EVALI(2) 
Z1=-EVALR(1)/(DSQRT(EVALI(1)**2+EVALR(1)**2)) 
Z2=-EVALR(2)/(DSQRT(EVALI(2)**2+EVALR(2)**2)) 
WN1=DSQRT(EVALI(1 )**2+EVALR(1 )**2) 
WN2=DSQRT(EVALI(2)**2+EVALR(2)**2) 
WRITE(37,25)Z1 ,Z2,WN1 ,WN2 
25   F0RMAT(4E15.5) 
500  CONTINUE 
RETURN 
END 
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B.l:   Effect  of  GM  on  Coupled  Model  Roll  Eigenvalues;  Real  Part 
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B.2:    Effect  of  GM  on  Coupled  Model  Roll  Eigenvalues;  Imaginary  Part 
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B.3:   Effect  of  GM  on  Coupled  Model  Upper  Steering  Eigenvalues 
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B.4:    Effect  of  GM  on  Coupled  Model  Lower  Steering  Eigenvalues 
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B.5:   Effect  of  Fn  on  Decoupled  Model  Damping  Ratio 
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GM  =  0.3  m  and   Fn  =  0.6 
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B.6:   Effect  of  Fn  on  Coupled  Model  Roll  Eigenvalues 
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B.7:    Effect  of  Fn  on  Coupled  Model  Upper  Steering  Eigenvalues 
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B.8:   Effect  of  Fn  on  Coupled  Model  Lower  Steering  Eigenvalues 
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B.9:     Effect  of  Fn  on  Coupled  Model  Roll  Eigenvalues 
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B.10:    Effect  of  Fn  on  Coupled  Model  Root  Locus 
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GM  =  0.3  m  and  Fn  =  0.3 
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B.ll:   Effect  of  Speed  Loss  on  Coupled  Model  Roll  Eigenvalues 
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B.12:   Effect  of  Speed  Loss  on  Coupled  Model  Root  Locus 
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B.13:     Effect  of  Speed  Loss  on  Coupled  Model  Upper  Steering 
Eigenvalues 
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B.14:     Effect  of  Speed  Loss  on  Coupled  Model  Lower  Steering 
Eigenvalues 
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B.15:   Effect  of  Speed  Loss  on  Coupled  Model  Sway  Velocity  with  GM  =  1.0  m 
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B.16:   Effect  of  Speed  Loss  on  Coupled  Model  Sway  Velocity  with  GM  =  3.0  m 
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B.17:   Effect  of  Speed  Loss  on  Coupled  Model  Roll  Angle  with  GM  =  1.0  m 
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B.18:   Effect  of  Speed  Loss  on  Coupled  Model  Roll  Angle  with  GM  =  3.0  m 
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B.19:   Effect  of  Speed  Loss  on  Coupled  Model  Root  Locus  with  GM  =  1.0  m 
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B.20:   Effect  of  Speed  Loss  on  Coupled  Model  Root  Locus  with  GM  =  3.0  m 
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